Click on tabs to display additional information.
# Load required libraries for statistical analysis and table creation
library(dplyr)
library(ggplot2)
library(gt)
library(broom)
library(car)
library(emmeans)
library(multcomp)
library(MASS) # For negative binomial regression
library(rstatix)
library(coin)
library(phyloseq)
library(microViz)
library(tidyverse)
# Define treatment order and color palette
treatment_order <- c(
"A- T- P-", # Control
"A- T- P+", # Parasite
"A+ T- P-", # Antibiotics
"A+ T- P+", # Antibiotics_Parasite
"A- T+ P-", # Temperature
"A- T+ P+", # Temperature_Parasite
"A+ T+ P-", # Antibiotics_Temperature
"A+ T+ P+" # Antibiotics_Temperature_Parasite
)
# Custom color palette matching treatment order
treatment_colors <- c(
"#1B9E77", # A- T- P- (Control)
"#D95F02", # A- T- P+ (Parasite)
"#7570B3", # A+ T- P- (Antibiotics)
"#E7298A", # A+ T- P+ (Antibiotics_Parasite)
"#66A61E", # A- T+ P- (Temperature)
"#E6AB02", # A- T+ P+ (Temperature_Parasite)
"#A6761D", # A+ T+ P- (Antibiotics_Temperature)
"#666666" # A+ T+ P+ (Antibiotics_Temperature_Parasite)
)
# Create named vector for color scale
treatment_color_scale <- setNames(treatment_colors, treatment_order)
# Function to extract sample data as dataframe from phyloseq object
samdatAsDataframe <- function(ps) {
samdat <- phyloseq::sample_data(ps)
df <- data.frame(samdat, check.names = FALSE, stringsAsFactors = FALSE)
return(df)
}
# Function to rename variables in phyloseq object
ps_rename <- function(ps, ...) {
ps <- microViz::ps_get(ps)
df <- samdatAsDataframe(ps)
df <- dplyr::rename(.data = df, ...)
phyloseq::sample_data(ps) <- df
return(ps)
}
# SourceFolder function
source(here::here("Code", "R", "Functions", "StartFunctions", "sourceFolder.R"))
# Import all helper functions found in `/Functions`
sourceFolder(here::here("Code", "R", "Functions", "StartFunctions"), T)
##
## Attaching package: 'scales'
## The following object is masked from 'package:purrr':
##
## discard
## The following object is masked from 'package:readr':
##
## col_factor
## The following object is masked from 'package:coin':
##
## pvalue
## Loading required package: future
##
## Attaching package: 'future'
## The following object is masked from 'package:survival':
##
## cluster
##
## Attaching package: 'ggpubr'
## The following object is masked from 'package:microViz':
##
## stat_chull
## Loading required package: grid
## ========================================
## ComplexHeatmap version 2.18.0
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
##
## If you use it in published research, please cite either one:
## - Gu, Z. Complex Heatmap Visualization. iMeta 2022.
## - Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional
## genomic data. Bioinformatics 2016.
##
##
## The new InteractiveComplexHeatmap package can directly export static
## complex heatmaps into an interactive Shiny app with zero effort. Have a try!
##
## This message can be suppressed by:
## suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================
## Loading required package: data.table
##
## Attaching package: 'data.table'
## The following objects are masked from 'package:lubridate':
##
## hour, isoweek, mday, minute, month, quarter, second, wday, week,
## yday, year
## The following object is masked from 'package:purrr':
##
## transpose
## The following objects are masked from 'package:dplyr':
##
## between, first, last
## Loading required package: doParallel
## Loading required package: foreach
##
## Attaching package: 'foreach'
## The following objects are masked from 'package:purrr':
##
## accumulate, when
## Loading required package: iterators
## Loading required package: GUniFrac
## Registered S3 method overwritten by 'rmutil':
## method from
## print.response httr
## Loading required package: magrittr
##
## Attaching package: 'magrittr'
## The following object is masked from 'package:purrr':
##
## set_names
## The following object is masked from 'package:tidyr':
##
## extract
## Loading required package: rlang
##
## Attaching package: 'rlang'
## The following object is masked from 'package:magrittr':
##
## set_names
## The following object is masked from 'package:data.table':
##
## :=
## The following objects are masked from 'package:purrr':
##
## %@%, flatten, flatten_chr, flatten_dbl, flatten_int, flatten_lgl,
## flatten_raw, invoke, splice
## Loading required package: vegan
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.6-4
##
## microbiome R package (microbiome.github.com)
##
##
##
## Copyright (C) 2011-2022 Leo Lahti,
## Sudarshan Shetty et al. <microbiome.github.io>
##
## Attaching package: 'microbiome'
## The following object is masked from 'package:vegan':
##
## diversity
## The following object is masked from 'package:scales':
##
## alpha
## The following object is masked from 'package:ggplot2':
##
## alpha
## The following object is masked from 'package:base':
##
## transform
## Loading required package: ape
##
## Attaching package: 'ape'
## The following object is masked from 'package:ggpubr':
##
## rotate
## The following object is masked from 'package:dplyr':
##
## where
## Loading required package: nlme
##
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
##
## collapse
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
##
## Attaching package: 'lme4'
## The following object is masked from 'package:nlme':
##
## lmList
## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## 3 files sourced from folder "/Users/michaelsieler/Dropbox/Mac (2)/Documents/Sharpton_Lab/Projects_Repository/Rules_of_Life/major-experiment-2023/Code/R/Functions/StartFunctions"
sourceFolder(here::here("Code", "R", "Functions", "HelperFunctions"), T)
## 9 files sourced from folder "/Users/michaelsieler/Dropbox/Mac (2)/Documents/Sharpton_Lab/Projects_Repository/Rules_of_Life/major-experiment-2023/Code/R/Functions/HelperFunctions"
# sourceFolder(here::here("Code", "R", "Functions", "AnalysisScripts"), T)
# Source the calculate_dynamic_metrics.R script to load required functions
source(here::here("Code", "R", "Functions", "AnalysisScripts", "calculate_dynamic_metrics.R"))
ps.tmp <- readRDS("/Users/michaelsieler/Dropbox/Mac (2)/Documents/Sharpton_Lab/Projects_Repository/Rules_of_Life/major-experiment-2023/Data/Robjects/pseq_uncleaned_05052025.rds")
ps.cleaned <-
ps.tmp %>%
## Update Metadata
ps_rename(Time = Timepoint) %>%
microViz::ps_mutate(
Treatment = case_when(
Antibiotics == 0 & Temperature == 0 & Pathogen == 0 ~ "A- T- P-",
Antibiotics == 0 & Temperature == 0 & Pathogen == 1 ~ "A- T- P+",
Antibiotics == 1 & Temperature == 0 & Pathogen == 0 ~ "A+ T- P-",
Antibiotics == 1 & Temperature == 0 & Pathogen == 1 ~ "A+ T- P+",
Antibiotics == 0 & Temperature == 1 & Pathogen == 0 ~ "A- T+ P-",
Antibiotics == 0 & Temperature == 1 & Pathogen == 1 ~ "A- T+ P+",
Antibiotics == 1 & Temperature == 1 & Pathogen == 0 ~ "A+ T+ P-",
Antibiotics == 1 & Temperature == 1 & Pathogen == 1 ~ "A+ T+ P+",
TRUE ~ "Unknown"
), .after = "Pathogen"
) %>%
microViz::ps_mutate(Sample = fecal.sample.number, .before = 1) %>%
microViz::ps_mutate(Sample = gsub("^f", "", Sample)) %>%
microViz::ps_filter(Treatment != "Unknown") %>%
microViz::ps_mutate(
History = case_when(
Antibiotics + Temperature == 0 ~ 0,
Antibiotics + Temperature == 1 ~ 1,
Antibiotics + Temperature == 2 ~ 2,
), .after = "Treatment"
) %>%
## Additional metadata updates, factorizing metadata
microViz::ps_mutate(
# Create treatment code
treatment_code = case_when(
Antibiotics == 0 & Temperature == 0 & Pathogen == 0 ~ "Aneg_Tneg_Pneg",
Antibiotics == 0 & Temperature == 0 & Pathogen == 1 ~ "Aneg_Tneg_Ppos",
Antibiotics == 1 & Temperature == 0 & Pathogen == 0 ~ "Apos_Tneg_Pneg",
Antibiotics == 1 & Temperature == 0 & Pathogen == 1 ~ "Apos_Tneg_Ppos",
Antibiotics == 0 & Temperature == 1 & Pathogen == 0 ~ "Aneg_Tpos_Pneg",
Antibiotics == 0 & Temperature == 1 & Pathogen == 1 ~ "Aneg_Tpos_Ppos",
Antibiotics == 1 & Temperature == 1 & Pathogen == 0 ~ "Apos_Tpos_Pneg",
Antibiotics == 1 & Temperature == 1 & Pathogen == 1 ~ "Apos_Tpos_Ppos"
),
# Create treatment group factor
treatment_group = case_when(
Antibiotics == 0 & Temperature == 0 & Pathogen == 1 ~ "Parasite",
Antibiotics == 1 & Temperature == 0 & Pathogen == 0 ~ "Antibiotics",
Antibiotics == 1 & Temperature == 0 & Pathogen == 1 ~ "Antibiotics_Parasite",
Antibiotics == 0 & Temperature == 1 & Pathogen == 0 ~ "Temperature",
Antibiotics == 0 & Temperature == 1 & Pathogen == 1 ~ "Temperature_Parasite",
Antibiotics == 1 & Temperature == 1 & Pathogen == 0 ~ "Antibiotics_Temperature",
Antibiotics == 1 & Temperature == 1 & Pathogen == 1 ~ "Antibiotics_Temperature_Parasite",
TRUE ~ "Control"
),
# Convert to factor with appropriate levels
treatment_group = factor(treatment_group,
levels = c("Control", "Parasite",
"Antibiotics", "Antibiotics_Parasite",
"Temperature", "Temperature_Parasite",
"Antibiotics_Temperature", "Antibiotics_Temperature_Parasite")
),
treatment_code = factor(treatment_code, levels = treatment_order),
# Create time point factor
time_point = factor(Time, levels = c(0, 14, 18, 25, 29, 60)),
# Create pathogen status factor
pathogen_status = factor(ifelse(Pathogen == 1, "Exposed", "Unexposed"),
levels = c("Unexposed", "Exposed")),
# Create sex factor
sex = factor(Sex, levels = c("M", "F"))
) %>%
microViz::ps_mutate(Treatment = factor(Treatment, levels = treatment_order)) %>%
microViz::ps_mutate(Exp_Type = case_when(
Treatment %in% c("A- T- P-", "A- T- P+") ~ "No prior stressor(s)",
Treatment %in% c("A+ T- P-", "A+ T- P+") ~ "Antibiotics",
Treatment %in% c("A- T+ P-", "A- T+ P+") ~ "Temperature",
Treatment %in% c("A+ T+ P-", "A+ T+ P+") ~ "Combined",
)) %>%
microViz::ps_mutate(Exp_Type = factor(Exp_Type, levels = c("No prior stressor(s)", "Antibiotics", "Temperature", "Combined"))) %>%
# Fix names for taxonomic ranks not identified
microViz::tax_fix(suffix_rank = "current", anon_unique = T, unknown = NA) %>%
# Filter for any samples that contain more than 5000 reads
microViz::ps_filter(sample_sums(.) > 5000) %>%
# Any taxa not found in at least 3 samples are removed
microViz::tax_filter(min_prevalence = 3, undetected = 0) %>%
# Remove any unwanted reads
microViz::tax_select(c("Mitochondria", "Chloroplast", "Eukaryota"), deselect = TRUE) %>%
microViz::tax_select(c("Bacteria, Phylum"), deselect = TRUE)
# ps.cleaned %>% microViz::samdat_tbl()
## Diversity Methods
diversity.method <- list()
diversity.method[["alpha"]] <- c("shannon",
"inverse_simpson",
"observed")
diversity.method[["beta"]] <- c("bray", "canberra")
## Stats/Plotting Variables
alpha.stats <- list() # save all stat results here
alpha.plots <- list() # save plots here
beta.stats <- list() # save all stat results here
beta.plots <- list() # save plots here
diffAbnd.stats <- list() # save all stat results here
diffAbnd.plots <- list() # save plots here
## Alpha
ps.cleaned_alpha <-
ps.cleaned %>% # phyloseq object we'll be using
microViz::ps_calc_diversity( # calculates various alpha diversity metrics
rank = "Genus", # What taxonomic rank do you want to calculate diversity metrics at
index = "shannon", # What diversity metric to use
varname = "Shannon__Genus", # Column name in sample data
exp = F # exponentiate the result or not
) %>%
microViz::ps_calc_diversity(
rank = "Genus",
index = "gini_simpson",
varname = "Simpson__Genus" # Column name in sample data
) %>%
microViz::ps_calc_richness( # related to ps_calc_diversity, this calculates richness or "observed" values
rank = "Genus",
varname = "Richness__Genus"
) %>%
# ps_calc_diversity.phy(
# varname = "Phylogenetic__Genus"
# ) %>% # Note: This is a helper function I made. You need to adjust function manually to change taxon rank (See: microViz_Helper.R)
# # ps_calc_diversity.phy() helper function does not properly name the column. Modify the following function as needed
# ps_rename("Phylogenetic__Genus" = "phylogenetic_Genus") %>% # Helper function (See: microViz_Helper.R)
ps_mutate(across(contains("__Genus"), norm_scores, .names = "{.col}_norm"))
##
## lambda W Shapiro.p.value A Anderson.p.value
## 420 0.475 0.987 3.787e-05 2.49 2.7e-06
##
## if (lambda > 0){TRANS = x ^ lambda}
## if (lambda == 0){TRANS = log(x)}
## if (lambda < 0){TRANS = -1 * x ^ lambda}
##
##
## lambda W Shapiro.p.value A Anderson.p.value
## 482 2.025 0.9761 2.845e-08 2.862 3.33e-07
##
## if (lambda > 0){TRANS = x ^ lambda}
## if (lambda == 0){TRANS = log(x)}
## if (lambda < 0){TRANS = -1 * x ^ lambda}
##
##
## lambda W Shapiro.p.value A Anderson.p.value
## 393 -0.2 0.9883 0.0001119 1.69 0.0002453
##
## if (lambda > 0){TRANS = x ^ lambda}
## if (lambda == 0){TRANS = log(x)}
## if (lambda < 0){TRANS = -1 * x ^ lambda}
ps.cleaned_alpha %>%
microViz::samdat_tbl()
## # A tibble: 594 × 37
## .sample_name Sample Tank.ID Time Antibiotics Temperature Pathogen Treatment
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <fct>
## 1 f1 1 1 0 0 0 0 A- T- P-
## 2 f2 2 1 0 0 0 0 A- T- P-
## 3 f3 3 1 0 0 0 0 A- T- P-
## 4 f31 31 2 0 0 0 0 A- T- P-
## 5 f32 32 2 0 0 0 0 A- T- P-
## 6 f33 33 2 0 0 0 0 A- T- P-
## 7 f61 61 3 0 0 0 0 A- T- P-
## 8 f62 62 3 0 0 0 0 A- T- P-
## 9 f63 63 3 0 0 0 0 A- T- P-
## 10 f91 91 4 0 0 0 1 A- T- P+
## # ℹ 584 more rows
## # ℹ 29 more variables: History <dbl>, Sex <chr>, Weight.g <dbl>,
## # Weight.mg <dbl>, Length.mm <dbl>, Body.Size <chr>, Mortality <dbl>,
## # Total.Worm.Count <dbl>, Flagged <dbl>, notes <chr>,
## # fecal.sample.number <chr>, gut.sample.number <chr>, adultworm.count <dbl>,
## # larvalworm.count <dbl>, notes.dissection <chr>, PLATE <chr>, WELL <chr>,
## # treatment_code <fct>, treatment_group <fct>, time_point <fct>, …
test.data.alpha <-
ps.cleaned_alpha %>%
microViz::samdat_tbl() %>%
dplyr::relocate(c(".sample_name"), .before = "Antibiotics") %>%
tidyr::pivot_longer(cols = contains("__Genus"), names_to = "Alpha.Metric", values_to = "Alpha.Score") %>%
dplyr::mutate(Data.Type = case_when(
Alpha.Metric %in% c("Shannon__Genus_norm", "Simpson__Genus_norm", "Richness__Genus_norm") ~ "Normalized",
TRUE ~ "Raw"
)) %>%
dplyr::relocate(c("Alpha.Metric", "Alpha.Score", "Data.Type"), .after = "Tank.ID") %>%
dplyr::relocate(c("Sample", "Treatment", "Time", "Tank.ID"), .before = 1)
# tidyr::pivot_longer(cols = contains("_norm"), names_to = "Alpha.Metric", values_to = "Alpha.Score") %>%
# dplyr::filter(Alpha.Metric %in% c("Simpson__Genus_norm", "Simpson__Genus_norm")) %>%
# tidyr::pivot_longer(cols = contains("__Genus"), names_to = "Raw.Metric", values_to = "Raw.Score")
#### Beta diversity matrices
##### All
beta.dist.mat <- list()
beta.dist.mat[["All"]] <- # saves the results of this loop to a list
furrr::future_map(diversity.method[["beta"]], function(x){ # Future_map runs the function in parallel
ps.cleaned_alpha %>% # is the phyloseq object we are looping over
microViz::tax_transform(trans = "identity", # tax transform transforms taxa counts
rank = ifelse(x == "gunifrac" || x == "wunifrac" || x == "unifrac", "unique", "Genus")) %>% # phylogenetic methods need "unique" ASV counts
microViz::dist_calc(x, gunifrac_alpha = 0.5) # calculates the distance between samples. gunifrac_alpha = 0.5 weights abundance into its calculation
}) %>%
setNames(diversity.method[["beta"]]) # assigns names of the beta methods to the list
beta.dist.mat[["All_60DPE"]] <- # saves the results of this loop to a list
furrr::future_map(diversity.method[["beta"]], function(x){ # Future_map runs the function in parallel
ps.cleaned_alpha %>% # is the phyloseq object we are looping over
microViz::ps_filter(Time == 60) %>%
microViz::tax_transform(trans = "identity", # tax transform transforms taxa counts
rank = ifelse(x == "gunifrac" || x == "wunifrac" || x == "unifrac", "unique", "Genus")) %>% # phylogenetic methods need "unique" ASV counts
microViz::dist_calc(x, gunifrac_alpha = 0.5) # calculates the distance between samples. gunifrac_alpha = 0.5 weights abundance into its calculation
}) %>%
setNames(diversity.method[["beta"]]) # assigns names of the beta methods to the list
##### Parasite exposure only
beta.dist.mat[["ParasiteExp"]] <- # saves the results of this loop to a list
furrr::future_map(diversity.method[["beta"]], function(x){ # Future_map runs the function in parallel
ps.cleaned_alpha %>% # is the phyloseq object we are looping over
microViz::ps_filter(Treatment %in% c(
"A- T- P-",
"A- T- P+"#,
# "A+ T- P-",
# "A+ T- P+",
# "A- T+ P-",
# "A- T+ P+",
# "A+ T+ P-",
# "A+ T+ P+"
)) %>%
microViz::tax_transform(trans = "identity", # tax transform transforms taxa counts
rank = ifelse(x == "gunifrac" || x == "wunifrac" || x == "unifrac", "unique", "Genus")) %>% # phylogenetic methods need "unique" ASV counts
microViz::dist_calc(x, gunifrac_alpha = 0.5) # calculates the distance between samples. gunifrac_alpha = 0.5 weights abundance into its calculation
}) %>%
setNames(diversity.method[["beta"]]) # assigns names of the beta methods to the list
beta.dist.mat[["ParasiteExp_60DPE"]] <- # saves the results of this loop to a list
furrr::future_map(diversity.method[["beta"]], function(x){ # Future_map runs the function in parallel
ps.cleaned_alpha %>% # is the phyloseq object we are looping over
microViz::ps_filter(Treatment %in% c(
"A- T- P-",
"A- T- P+"#,
# "A+ T- P-",
# "A+ T- P+",
# "A- T+ P-",
# "A- T+ P+",
# "A+ T+ P-",
# "A+ T+ P+"
)) %>%
microViz::ps_filter(Time == 60) %>%
microViz::tax_transform(trans = "identity", # tax transform transforms taxa counts
rank = ifelse(x == "gunifrac" || x == "wunifrac" || x == "unifrac", "unique", "Genus")) %>% # phylogenetic methods need "unique" ASV counts
microViz::dist_calc(x, gunifrac_alpha = 0.5) # calculates the distance between samples. gunifrac_alpha = 0.5 weights abundance into its calculation
}) %>%
setNames(diversity.method[["beta"]]) # assigns names of the beta methods to the list
##### Prior stressors and parasite exposure
beta.dist.mat[["PriorStressParaExp"]] <- # saves the results of this loop to a list
furrr::future_map(diversity.method[["beta"]], function(x){ # Future_map runs the function in parallel
ps.cleaned_alpha %>% # is the phyloseq object we are looping over
microViz::ps_filter(Treatment %in% c(
# "A- T- P-",
"A- T- P+",
# "A+ T- P-",
"A+ T- P+",
# "A- T+ P-",
"A- T+ P+",
# "A+ T+ P-",
"A+ T+ P+"
)) %>%
microViz::tax_transform(trans = "identity", # tax transform transforms taxa counts
rank = ifelse(x == "gunifrac" || x == "wunifrac" || x == "unifrac", "unique", "Genus")) %>% # phylogenetic methods need "unique" ASV counts
microViz::dist_calc(x, gunifrac_alpha = 0.5) # calculates the distance between samples. gunifrac_alpha = 0.5 weights abundance into its calculation
}) %>%
setNames(diversity.method[["beta"]]) # assigns names of the beta methods to the list
beta.dist.mat[["PriorStressParaExp_60DPE"]] <- # saves the results of this loop to a list
furrr::future_map(diversity.method[["beta"]], function(x){ # Future_map runs the function in parallel
ps.cleaned_alpha %>% # is the phyloseq object we are looping over
microViz::ps_filter(Treatment %in% c(
# "A- T- P-",
"A- T- P+",
# "A+ T- P-",
"A+ T- P+",
# "A- T+ P-",
"A- T+ P+",
# "A+ T+ P-",
"A+ T+ P+"
)) %>%
microViz::ps_filter(Time == 60) %>%
microViz::tax_transform(trans = "identity", # tax transform transforms taxa counts
rank = ifelse(x == "gunifrac" || x == "wunifrac" || x == "unifrac", "unique", "Genus")) %>% # phylogenetic methods need "unique" ASV counts
microViz::dist_calc(x, gunifrac_alpha = 0.5) # calculates the distance between samples. gunifrac_alpha = 0.5 weights abundance into its calculation
}) %>%
setNames(diversity.method[["beta"]]) # assigns names of the beta methods to the list
##### Prior stressors and NO parasite exposure
beta.dist.mat[["PriorStressNoParaExp"]] <- # saves the results of this loop to a list
furrr::future_map(diversity.method[["beta"]], function(x){ # Future_map runs the function in parallel
ps.cleaned_alpha %>% # is the phyloseq object we are looping over
microViz::ps_filter(Treatment %in% c(
"A- T- P-",
# "A- T- P+",
"A+ T- P-",
# "A+ T- P+",
"A- T+ P-",
# "A- T+ P+",
"A+ T+ P-"#,
# "A+ T+ P+"
)) %>%
microViz::tax_transform(trans = "identity", # tax transform transforms taxa counts
rank = ifelse(x == "gunifrac" || x == "wunifrac" || x == "unifrac", "unique", "Genus")) %>% # phylogenetic methods need "unique" ASV counts
microViz::dist_calc(x, gunifrac_alpha = 0.5) # calculates the distance between samples. gunifrac_alpha = 0.5 weights abundance into its calculation
}) %>%
setNames(diversity.method[["beta"]]) # assigns names of the beta methods to the list
beta.dist.mat[["PriorStressNoParaExp_60DPE"]] <- # saves the results of this loop to a list
furrr::future_map(diversity.method[["beta"]], function(x){ # Future_map runs the function in parallel
ps.cleaned_alpha %>% # is the phyloseq object we are looping over
microViz::ps_filter(Treatment %in% c(
"A- T- P-",
# "A- T- P+",
"A+ T- P-",
# "A+ T- P+",
"A- T+ P-",
# "A- T+ P+",
"A+ T+ P-"#,
# "A+ T+ P+"
)) %>%
microViz::ps_filter(Time == 60) %>%
microViz::tax_transform(trans = "identity", # tax transform transforms taxa counts
rank = ifelse(x == "gunifrac" || x == "wunifrac" || x == "unifrac", "unique", "Genus")) %>% # phylogenetic methods need "unique" ASV counts
microViz::dist_calc(x, gunifrac_alpha = 0.5) # calculates the distance between samples. gunifrac_alpha = 0.5 weights abundance into its calculation
}) %>%
setNames(diversity.method[["beta"]]) # assigns names of the beta methods to the list
#### Beta Personalization Scores
data.beta.list <- list()
# Save dataframe
data.beta.list[["All"]] <-
ps.cleaned_alpha %>%
# Calculate Distance matrix
microViz::dist_calc("bray") %>%
# Get distance matric
microViz::dist_get() %>%
# Pat Schloss methods to view beta diversity matrix data
# Source: https://youtu.be/Gm-kg3TuML4?si=nYLLAxm7uv8aYG27&t=250
# Convert into a matrix object
as.matrix() %>%
# Convert matrix into a tibble
tibble::as_tibble(rownames = "Sample_a") %>%
# Convert dataframe into "long" format
tidyr::pivot_longer(-Sample_a, names_to = "Sample_b", values_to = "Dist") %>%
# Remove any rows where names equal eachother
dplyr::filter(Sample_a != Sample_b) %>%
# Merge sample dataframe
dplyr::right_join((ps.cleaned_alpha %>%
microViz::samdat_tbl() %>% ungroup() %>% arrange(Time, Treatment, Tank.ID) %>%
dplyr::select(!contains("_Genus")) %>%
dplyr::select(fecal.sample.number, Time, Treatment, Tank.ID)
),
by = c("Sample_a" = "fecal.sample.number")) %>%
dplyr::right_join((ps.cleaned_alpha %>%
microViz::samdat_tbl() %>%
dplyr::ungroup() %>%
dplyr::arrange(Time, Treatment, Tank.ID) %>%
dplyr::select(!contains("_Genus")) %>%
dplyr::select(fecal.sample.number, Time, Treatment, Tank.ID)
),
by = c("Sample_b" = "fecal.sample.number"), suffix = c("_a", "_b")) %>% # add suffix to differentiate between sample metadata
# Move columns around
dplyr::select(order(colnames(.))) %>%
dplyr::relocate(dplyr::any_of(c("Sample_a", "Sample_b", "Dist")), .before = 1) %>%
dplyr::relocate(dplyr::contains(c("Tank")), .after = dplyr::last_col()) %>%
dplyr::ungroup() %>%
# Calculate the median distance of samples to all other samples (Measure of "personalization")
dplyr::group_by(Sample_a, Treatment_a, Time_a, Tank.ID_a) %>%
dplyr::summarize(Beta.Score = median(Dist), #*100,
Beta.Score.Mean = mean(Dist), #*100,
count = dplyr::n()) %>%
dplyr::ungroup() %>%
# Rename Columns
dplyr::rename(fecal.sample.number = Sample_a,
Treatment = Treatment_a,
Time = Time_a,
Tank.ID = Tank.ID_a) %>%
dplyr::right_join(ps.cleaned_alpha %>%
microViz::samdat_tbl() %>%
dplyr::select(!contains("_Genus")),
by = c("fecal.sample.number", "Time", "Treatment","Tank.ID")) %>%
dplyr::ungroup() %>%
# Normalize diversity scores
dplyr::mutate(dplyr::across(dplyr::contains("Beta.Score"), norm_scores, .names = "{.col}_norm"), .after = "Beta.Score") %>%
# Create a column to save what type of beta metric was calculated
dplyr::mutate(Beta.Metric = "Bray-Curtis", .before = "Beta.Score") %>%
dplyr::mutate(Sample = gsub("^f", "", fecal.sample.number)) %>%
dplyr::relocate(Sample, .before = 1) %>%
dplyr::relocate(fecal.sample.number, .after = 36)
## `summarise()` has grouped output by 'Sample_a', 'Treatment_a', 'Time_a'. You
## can override using the `.groups` argument.
##
## lambda W Shapiro.p.value A Anderson.p.value
## 291 -2.75 0.9882 0.0001004 1.05 0.009232
##
## if (lambda > 0){TRANS = x ^ lambda}
## if (lambda == 0){TRANS = log(x)}
## if (lambda < 0){TRANS = -1 * x ^ lambda}
##
##
## lambda W Shapiro.p.value A Anderson.p.value
## 235 -4.15 0.9883 0.0001067 1.004 0.01196
##
## if (lambda > 0){TRANS = x ^ lambda}
## if (lambda == 0){TRANS = log(x)}
## if (lambda < 0){TRANS = -1 * x ^ lambda}
data.beta.list[["All"]] <- data.beta.list[["All"]] %>%
dplyr::full_join(.,
ps.cleaned_alpha %>%
# tax_agg("Genus") %>%
microViz::dist_calc("canberra") %>%
# Get distance matric
microViz::dist_get() %>%
# Pat Schloss methods to view beta diversity matrix data
# Source: https://youtu.be/Gm-kg3TuML4?si=nYLLAxm7uv8aYG27&t=250
# Convert into a matrix object
as.matrix() %>%
# Convert matrix into a tibble
tibble::as_tibble(rownames = "Sample_a") %>%
# Convert dataframe into "long" format
tidyr::pivot_longer(-Sample_a, names_to = "Sample_b", values_to = "Dist") %>%
# Remove any rows where names equal eachother
dplyr::filter(Sample_a != Sample_b) %>%
# Merge sample dataframe
dplyr::right_join((ps.cleaned_alpha %>%
microViz::samdat_tbl() %>% ungroup() %>% arrange(Time, Treatment, Tank.ID) %>%
dplyr::select(!contains("_Genus")) %>%
dplyr::select(fecal.sample.number, Time, Treatment, Tank.ID)
),
by = c("Sample_a" = "fecal.sample.number")) %>%
dplyr::right_join((ps.cleaned_alpha %>%
microViz::samdat_tbl() %>%
dplyr::ungroup() %>%
dplyr::arrange(Time, Treatment, Tank.ID) %>%
dplyr::select(!contains("_Genus")) %>%
dplyr::select(fecal.sample.number, Time, Treatment, Tank.ID)
),
by = c("Sample_b" = "fecal.sample.number"), suffix = c("_a", "_b")) %>% # add suffix to differentiate between sample metadata
# Move columns around
dplyr::select(order(colnames(.))) %>%
dplyr::relocate(dplyr::any_of(c("Sample_a", "Sample_b", "Dist")), .before = 1) %>%
dplyr::relocate(dplyr::contains(c("Tank")), .after = dplyr::last_col()) %>%
dplyr::ungroup() %>%
# Calculate the median distance of samples to all other samples (Measure of "personalization")
dplyr::group_by(Sample_a, Treatment_a, Time_a, Tank.ID_a) %>%
dplyr::summarize(Beta.Score = median(Dist), #*100,
Beta.Score.Mean = mean(Dist), #*100,
count = dplyr::n()) %>%
dplyr::ungroup() %>%
# Rename Columns
dplyr::rename(fecal.sample.number = Sample_a,
Treatment = Treatment_a,
Time = Time_a,
Tank.ID = Tank.ID_a) %>%
dplyr::right_join(ps.cleaned_alpha %>%
microViz::samdat_tbl() %>%
dplyr::select(!contains("_Genus")),
by = c("fecal.sample.number", "Time", "Treatment","Tank.ID")) %>%
dplyr::ungroup() %>%
# Normalize diversity scores
dplyr::mutate(dplyr::across(contains("Beta.Score"), norm_scores, .names = "{.col}_norm"), .after = "Beta.Score") %>%
# Create a column to save what type of beta metric was calculated
dplyr::mutate(Beta.Metric = "Canberra", .before = "Beta.Score")
)
## `summarise()` has grouped output by 'Sample_a', 'Treatment_a', 'Time_a'. You
## can override using the `.groups` argument.
##
## lambda W Shapiro.p.value A Anderson.p.value
## 76 -8.125 0.9925 0.004258 0.6477 0.09061
##
## if (lambda > 0){TRANS = x ^ lambda}
## if (lambda == 0){TRANS = log(x)}
## if (lambda < 0){TRANS = -1 * x ^ lambda}
##
##
## lambda W Shapiro.p.value A Anderson.p.value
## 51 -8.75 0.9948 0.04029 0.4685 0.2483
##
## if (lambda > 0){TRANS = x ^ lambda}
## if (lambda == 0){TRANS = log(x)}
## if (lambda < 0){TRANS = -1 * x ^ lambda}
## Joining with `by = join_by(Sample, Treatment, Time, Tank.ID, Beta.Metric,
## Beta.Score, Beta.Score_norm, Beta.Score.Mean_norm, Beta.Score.Mean, count,
## .sample_name, Antibiotics, Temperature, Pathogen, History, Sex, Weight.g,
## Weight.mg, Length.mm, Body.Size, Mortality, Total.Worm.Count, Flagged, notes,
## gut.sample.number, adultworm.count, larvalworm.count, notes.dissection, PLATE,
## WELL, treatment_code, treatment_group, time_point, pathogen_status, sex,
## fecal.sample.number, Exp_Type)`
# Prepare beta diversity data
test.data.beta <-
data.beta.list[["All"]] %>%
tidyr::pivot_longer(cols = c("Beta.Score", "Beta.Score_norm"), names_to = "Beta.Metric.Type", values_to = "Beta.Score") %>%
dplyr::relocate(c("Beta.Metric.Type", "Beta.Score"), .after = "Beta.Metric") %>%
dplyr::mutate(Data.Type = case_when(
Beta.Metric.Type %in% c(" Beta.Score", "Beta.Score_norm") ~ "Normalized",
TRUE ~ "Raw"
), .after = "Beta.Score") %>%
# Remove the temporary Beta.Metric.Type column
dplyr::select(-c(Beta.Metric.Type, count)) %>%
dplyr::mutate(Beta.Metric = ifelse(Data.Type == "Normalized", paste0(Beta.Metric, "_norm"), Beta.Metric)) %>%
dplyr::arrange(Sample, Time, Tank.ID, Data.Type)
# Calculate all metrics at once using the combined function
test.alpha.all_metrics <- calculate_all_dynamic_metrics(
data = test.data.alpha,
metric_col = "Alpha.Metric",
score_col = "Alpha.Score",
ref_group = NULL, # Using all groups as reference
time_col = "Time",
group_col = "Treatment",
use_mean = TRUE
)
## Joining with `by = join_by(Sample, Time, Treatment, Tank.ID, Alpha.Metric,
## Data.Type, Displacement)`
## Joining with `by = join_by(Sample, Time, Treatment, Tank.ID, Alpha.Metric,
## Data.Type)`
# Calculate all metrics at once using the combined function
test.beta.all_metrics <- calculate_all_dynamic_metrics(
data = test.data.beta,
metric_col = "Beta.Metric",
score_col = "Beta.Score",
ref_group = NULL,
time_col = "Time",
group_col = "Treatment",
use_mean = TRUE
)
## Joining with `by = join_by(Sample, Time, Treatment, Tank.ID, Beta.Metric,
## Data.Type, Displacement)`
## Joining with `by = join_by(Sample, Time, Treatment, Tank.ID, Beta.Metric,
## Data.Type)`
# Merge Displacement data with Alpha.Metric and Alpha.Score
data.alpha.merged <-
test.alpha.all_metrics$volatility %>%
# dplyr::filter()
# First pivot the SADMM metrics
# dplyr::group_by(Data.Type) %>%
tidyr::pivot_longer(
cols = c("Alpha.Score", "Ref.Score", "Displacement", "Velocity", "Acceleration", "Volatility"),
names_to = "METRIC",
values_to = "SCORE"
) %>%
dplyr::mutate(METRIC = factor(METRIC, levels = c("Alpha.Score", "Ref.Score", "Displacement", "Velocity", "Acceleration", "Volatility")),
Data.Type = factor(Data.Type, levels = c("Raw", "Normalized")))
# Merge Displacement data with Alpha.Metric and Alpha.Score
data.beta.merged <-
test.beta.all_metrics$volatility %>%
# dplyr::filter()
# First pivot the SADMM metrics
# dplyr::group_by(Data.Type) %>%
tidyr::pivot_longer(
cols = c("Beta.Score", "Ref.Score", "Displacement", "Velocity", "Acceleration", "Volatility"),
names_to = "METRIC",
values_to = "SCORE"
) %>%
dplyr::mutate(METRIC = factor(METRIC, levels = c("Beta.Score", "Ref.Score", "Displacement", "Velocity", "Acceleration", "Volatility")),
Data.Type = factor(Data.Type, levels = c("Raw", "Normalized")))
# Set the subsection for this analysis
tmp.resSubSection <- "All"
# GLM Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]] <-
test.alpha.all_metrics$volatility %>%
# Clean up cell value names by removing any strings including and after "__"
cutCellNames(col = "Alpha.Metric", sep = "__") %>%
# Filter out raw data, leave only normalized data
dplyr::filter(Data.Type != "Raw") %>%
# Run GLM models
run_glm_models(formula_str = "Displacement ~ Treatment",
alpha_metric_col = "Alpha.Metric",
alpha_score_col = "Displacement",
family_str = "gaussian")
# ANOVA Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["ANOVA"]] <-
run_glm_anova(alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]])
# Tukey Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["Tukey"]] <-
test.alpha.all_metrics$volatility %>%
# Clean up cell value names by removing any strings including and after "__"
cutCellNames(col = "Alpha.Metric", sep = "__") %>%
# Filter out raw data, leave only normalized data
dplyr::filter(Data.Type != "Raw") %>%
# Run Tukey test on GLM models
run_tukey_glm(., "Displacement", "Alpha.Metric", c("Treatment"),
family_str = "gaussian",
group_by_var = NULL)
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
| GLM Results | |||||
| glm(Displacement ~ Treatment); All treatments | |||||
| term | estimate | std.error | statistic | p.value | p.adj.sig |
|---|---|---|---|---|---|
| Shannon | |||||
| (Intercept) | 0.078 | 0.019 | 3.995 | <0.001 | **** |
| A- T- P+ | −0.045 | 0.027 | −1.627 | 0.104 | ns |
| A+ T- P- | −0.051 | 0.028 | −1.815 | 0.070 | ns |
| A+ T- P+ | −0.071 | 0.028 | −2.584 | 0.010 | ** |
| A- T+ P- | −0.004 | 0.028 | −0.151 | ≥0.25 | ns |
| A- T+ P+ | −0.092 | 0.028 | −3.302 | 0.001 | ** |
| A+ T+ P- | −0.113 | 0.030 | −3.799 | <0.001 | *** |
| A+ T+ P+ | −0.061 | 0.029 | −2.149 | 0.032 | * |
| Simpson | |||||
| (Intercept) | 0.042 | 0.028 | 1.527 | 0.127 | ns |
| A- T- P+ | −0.049 | 0.039 | −1.238 | 0.216 | ns |
| A+ T- P- | −0.034 | 0.040 | −0.841 | ≥0.25 | ns |
| A+ T- P+ | −0.085 | 0.039 | −2.165 | 0.031 | * |
| A- T+ P- | 0.035 | 0.041 | 0.874 | ≥0.25 | ns |
| A- T+ P+ | −0.154 | 0.040 | −3.867 | <0.001 | *** |
| A+ T+ P- | −0.181 | 0.042 | −4.270 | <0.001 | **** |
| A+ T+ P+ | −0.091 | 0.041 | −2.232 | 0.026 | * |
| Richness | |||||
| (Intercept) | 0.221 | 0.024 | 9.249 | <0.001 | **** |
| A- T- P+ | 0.009 | 0.034 | 0.265 | ≥0.25 | ns |
| A+ T- P- | −0.049 | 0.035 | −1.430 | 0.153 | ns |
| A+ T- P+ | −0.045 | 0.034 | −1.325 | 0.186 | ns |
| A- T+ P- | −0.016 | 0.035 | −0.467 | ≥0.25 | ns |
| A- T+ P+ | 0.016 | 0.034 | 0.456 | ≥0.25 | ns |
| A+ T+ P- | 0.029 | 0.037 | 0.783 | ≥0.25 | ns |
| A+ T+ P+ | −0.054 | 0.035 | −1.533 | 0.126 | ns |
| ANOVA of GLM | ||||
| ANOVA(GLM(Displacement ~ Treatment), type = 2); All treatments | ||||
| term | statistic | df | p.value | sig |
|---|---|---|---|---|
| Shannon | ||||
| Treatment | 25.086 | 7.000 | <0.001 | *** |
| Simpson | ||||
| Treatment | 42.565 | 7.000 | <0.001 | **** |
| Richness | ||||
| Treatment | 11.161 | 7.000 | 0.132 | ns |
| Pairwise Tukey's HSD, p.adj: Dunnett | ||||||||
| Tukey(Displacement ~ Treatment); All treatments | ||||||||
| term | .y. | group1 | group2 | estimate | std.error | statistic | adj.p.value | Variable |
|---|---|---|---|---|---|---|---|---|
| Shannon | ||||||||
| Treatment | Alpha.Score | A- T- P+ | A- T- P- | −0.045 | 0.027 | −1.627 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T- P- | A- T- P- | −0.051 | 0.028 | −1.815 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T- P+ | A- T- P- | −0.071 | 0.028 | −2.584 | 0.161 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A- T- P- | −0.004 | 0.028 | −0.151 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A- T- P- | −0.092 | 0.028 | −3.302 | 0.022 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T- P- | −0.113 | 0.030 | −3.799 | 0.004 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T- P- | −0.061 | 0.029 | −2.149 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T- P- | A- T- P+ | −0.006 | 0.028 | −0.220 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T- P+ | A- T- P+ | −0.027 | 0.028 | −0.968 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A- T- P+ | 0.040 | 0.028 | 1.421 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A- T- P+ | −0.047 | 0.028 | −1.695 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T- P+ | −0.068 | 0.030 | −2.298 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T- P+ | −0.017 | 0.029 | −0.582 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T- P+ | A+ T- P- | −0.021 | 0.028 | −0.730 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A+ T- P- | 0.047 | 0.029 | 1.608 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A+ T- P- | −0.041 | 0.028 | −1.446 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A+ T- P- | −0.062 | 0.030 | −2.056 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A+ T- P- | −0.010 | 0.029 | −0.360 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A+ T- P+ | 0.067 | 0.029 | 2.348 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A+ T- P+ | −0.020 | 0.028 | −0.729 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A+ T- P+ | −0.042 | 0.030 | −1.392 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A+ T- P+ | 0.010 | 0.029 | 0.353 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A- T+ P- | −0.088 | 0.029 | −3.044 | 0.048 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T+ P- | −0.109 | 0.031 | −3.549 | 0.009 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T+ P- | −0.057 | 0.029 | −1.936 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T+ P+ | −0.021 | 0.030 | −0.706 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T+ P+ | 0.031 | 0.029 | 1.058 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A+ T+ P- | 0.052 | 0.031 | 1.684 | ≥0.25 | Treatment |
| Simpson | ||||||||
| Treatment | Alpha.Score | A- T- P+ | A- T- P- | −0.049 | 0.039 | −1.238 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T- P- | A- T- P- | −0.034 | 0.040 | −0.841 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T- P+ | A- T- P- | −0.085 | 0.039 | −2.165 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A- T- P- | 0.035 | 0.041 | 0.874 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A- T- P- | −0.154 | 0.040 | −3.867 | 0.003 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T- P- | −0.181 | 0.042 | −4.270 | <0.001 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T- P- | −0.091 | 0.041 | −2.232 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T- P- | A- T- P+ | 0.015 | 0.040 | 0.372 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T- P+ | A- T- P+ | −0.037 | 0.039 | −0.935 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A- T- P+ | 0.084 | 0.041 | 2.070 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A- T- P+ | −0.105 | 0.040 | −2.645 | 0.140 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T- P+ | −0.133 | 0.042 | −3.128 | 0.038 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T- P+ | −0.042 | 0.041 | −1.041 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T- P+ | A+ T- P- | −0.052 | 0.040 | −1.287 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A+ T- P- | 0.069 | 0.041 | 1.672 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A+ T- P- | −0.120 | 0.040 | −2.963 | 0.060 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A+ T- P- | −0.148 | 0.043 | −3.421 | 0.015 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A+ T- P- | −0.057 | 0.041 | −1.381 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A+ T- P+ | 0.121 | 0.041 | 2.962 | 0.061 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A+ T- P+ | −0.068 | 0.040 | −1.706 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A+ T- P+ | −0.096 | 0.043 | −2.248 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A+ T- P+ | −0.005 | 0.041 | −0.134 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A- T+ P- | −0.189 | 0.041 | −4.604 | <0.001 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T+ P- | −0.217 | 0.044 | −4.957 | <0.001 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T+ P- | −0.126 | 0.042 | −3.006 | 0.053 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T+ P+ | −0.028 | 0.043 | −0.649 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T+ P+ | 0.063 | 0.041 | 1.521 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A+ T+ P- | 0.091 | 0.044 | 2.063 | ≥0.25 | Treatment |
| Richness | ||||||||
| Treatment | Alpha.Score | A- T- P+ | A- T- P- | 0.009 | 0.034 | 0.265 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T- P- | A- T- P- | −0.049 | 0.035 | −1.430 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T- P+ | A- T- P- | −0.045 | 0.034 | −1.325 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A- T- P- | −0.016 | 0.035 | −0.467 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A- T- P- | 0.016 | 0.034 | 0.456 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T- P- | 0.029 | 0.037 | 0.783 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T- P- | −0.054 | 0.035 | −1.533 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T- P- | A- T- P+ | −0.058 | 0.035 | −1.691 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T- P+ | A- T- P+ | −0.054 | 0.034 | −1.589 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A- T- P+ | −0.025 | 0.035 | −0.724 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A- T- P+ | 0.007 | 0.034 | 0.194 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T- P+ | 0.020 | 0.037 | 0.538 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T- P+ | −0.063 | 0.035 | −1.789 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T- P+ | A+ T- P- | 0.004 | 0.035 | 0.122 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A+ T- P- | 0.033 | 0.036 | 0.925 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A+ T- P- | 0.065 | 0.035 | 1.860 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A+ T- P- | 0.078 | 0.037 | 2.093 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A+ T- P- | −0.005 | 0.036 | −0.126 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A+ T- P+ | 0.029 | 0.035 | 0.816 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A+ T- P+ | 0.061 | 0.035 | 1.761 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A+ T- P+ | 0.074 | 0.037 | 2.002 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A+ T- P+ | −0.009 | 0.035 | −0.248 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A- T+ P- | 0.032 | 0.035 | 0.903 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T+ P- | 0.045 | 0.038 | 1.193 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T+ P- | −0.038 | 0.036 | −1.034 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T+ P+ | 0.013 | 0.037 | 0.353 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T+ P+ | −0.070 | 0.036 | −1.955 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A+ T+ P- | −0.083 | 0.038 | −2.180 | ≥0.25 | Treatment |
# Set the subsection for this analysis
tmp.resSubSection <- "All_60DPE"
# GLM Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]] <-
test.alpha.all_metrics$displacement %>%
dplyr::filter(Time == 60) %>%
# Clean up cell value names by removing any strings including and after "__"
cutCellNames(col = "Alpha.Metric", sep = "__") %>%
# Filter out raw data, leave only normalized data
dplyr::filter(Data.Type != "Raw") %>%
# Run GLM models
run_glm_models(formula_str = "Displacement ~ Treatment",
alpha_metric_col = "Alpha.Metric",
alpha_score_col = "Displacement",
family_str = "gaussian")
# ANOVA Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["ANOVA"]] <-
run_glm_anova(alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]])
# Tukey Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["Tukey"]] <-
test.alpha.all_metrics$displacement %>%
dplyr::filter(Time == 60) %>%
# Clean up cell value names by removing any strings including and after "__"
cutCellNames(col = "Alpha.Metric", sep = "__") %>%
# Filter out raw data, leave only normalized data
dplyr::filter(Data.Type != "Raw") %>%
# Run Tukey test on GLM models
run_tukey_glm(., "Displacement", "Alpha.Metric", c("Treatment"),
family_str = "gaussian",
group_by_var = NULL)
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
| GLM Results - 60 DPE | |||||
| glm(Displacement ~ Treatment); All treatments at 60 DPE | |||||
| term | estimate | std.error | statistic | p.value | p.adj.sig |
|---|---|---|---|---|---|
| Shannon | |||||
| (Intercept) | 0.091 | 0.037 | 2.457 | 0.015 | * |
| A- T- P+ | −0.082 | 0.052 | −1.590 | 0.113 | ns |
| A+ T- P- | −0.034 | 0.055 | −0.617 | ≥0.25 | ns |
| A+ T- P+ | −0.083 | 0.052 | −1.590 | 0.113 | ns |
| A- T+ P- | −0.041 | 0.057 | −0.719 | ≥0.25 | ns |
| A- T+ P+ | −0.159 | 0.054 | −2.965 | 0.003 | ** |
| A+ T+ P- | −0.191 | 0.068 | −2.815 | 0.005 | ** |
| A+ T+ P+ | −0.079 | 0.058 | −1.364 | 0.174 | ns |
| Simpson | |||||
| (Intercept) | 0.008 | 0.049 | 0.162 | ≥0.25 | ns |
| A- T- P+ | −0.070 | 0.069 | −1.011 | ≥0.25 | ns |
| A+ T- P- | 0.023 | 0.073 | 0.311 | ≥0.25 | ns |
| A+ T- P+ | −0.070 | 0.070 | −1.002 | ≥0.25 | ns |
| A- T+ P- | 0.030 | 0.076 | 0.399 | ≥0.25 | ns |
| A- T+ P+ | −0.249 | 0.072 | −3.483 | <0.001 | *** |
| A+ T+ P- | −0.332 | 0.091 | −3.665 | <0.001 | *** |
| A+ T+ P+ | −0.106 | 0.077 | −1.387 | 0.167 | ns |
| Richness | |||||
| (Intercept) | 0.338 | 0.032 | 10.566 | <0.001 | **** |
| A- T- P+ | 0.025 | 0.045 | 0.552 | ≥0.25 | ns |
| A+ T- P- | −0.064 | 0.047 | −1.349 | 0.179 | ns |
| A+ T- P+ | −0.108 | 0.046 | −2.363 | 0.019 | * |
| A- T+ P- | −0.089 | 0.049 | −1.804 | 0.073 | ns |
| A- T+ P+ | 0.022 | 0.047 | 0.462 | ≥0.25 | ns |
| A+ T+ P- | 0.067 | 0.059 | 1.141 | ≥0.25 | ns |
| A+ T+ P+ | −0.028 | 0.050 | −0.560 | ≥0.25 | ns |
| ANOVA of GLM - 60 DPE | ||||
| ANOVA(GLM(Displacement ~ Treatment), type = 2); All treatments at 60 DPE | ||||
| term | statistic | df | p.value | sig |
|---|---|---|---|---|
| Shannon | ||||
| Treatment | 14.694 | 7.000 | 0.040 | * |
| Simpson | ||||
| Treatment | 31.743 | 7.000 | <0.001 | **** |
| Richness | ||||
| Treatment | 19.090 | 7.000 | 0.008 | ** |
| Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE | ||||||||
| Tukey(Displacement ~ Treatment); All treatments at 60 DPE | ||||||||
| term | .y. | group1 | group2 | estimate | std.error | statistic | adj.p.value | Variable |
|---|---|---|---|---|---|---|---|---|
| Shannon | ||||||||
| Treatment | Alpha.Score | A- T- P+ | A- T- P- | −0.082 | 0.052 | −1.590 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T- P- | A- T- P- | −0.034 | 0.055 | −0.617 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T- P+ | A- T- P- | −0.083 | 0.052 | −1.590 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A- T- P- | −0.041 | 0.057 | −0.719 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A- T- P- | −0.159 | 0.054 | −2.965 | 0.059 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T- P- | −0.191 | 0.068 | −2.815 | 0.089 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T- P- | −0.079 | 0.058 | −1.364 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T- P- | A- T- P+ | 0.049 | 0.054 | 0.894 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T- P+ | A- T- P+ | −0.001 | 0.052 | −0.022 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A- T- P+ | 0.041 | 0.057 | 0.731 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A- T- P+ | −0.077 | 0.053 | −1.442 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T- P+ | −0.109 | 0.068 | −1.611 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T- P+ | 0.004 | 0.057 | 0.066 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T- P+ | A+ T- P- | −0.050 | 0.055 | −0.904 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A+ T- P- | −0.007 | 0.059 | −0.122 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A+ T- P- | −0.126 | 0.056 | −2.234 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A+ T- P- | −0.158 | 0.070 | −2.254 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A+ T- P- | −0.045 | 0.060 | −0.748 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A+ T- P+ | 0.043 | 0.057 | 0.743 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A+ T- P+ | −0.076 | 0.054 | −1.402 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A+ T- P+ | −0.108 | 0.068 | −1.581 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A+ T- P+ | 0.005 | 0.058 | 0.085 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A- T+ P- | −0.118 | 0.058 | −2.027 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T+ P- | −0.150 | 0.072 | −2.097 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T+ P- | −0.038 | 0.062 | −0.607 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T+ P+ | −0.032 | 0.069 | −0.463 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T+ P+ | 0.081 | 0.059 | 1.368 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A+ T+ P- | 0.113 | 0.072 | 1.562 | ≥0.25 | Treatment |
| Simpson | ||||||||
| Treatment | Alpha.Score | A- T- P+ | A- T- P- | −0.070 | 0.069 | −1.011 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T- P- | A- T- P- | 0.023 | 0.073 | 0.311 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T- P+ | A- T- P- | −0.070 | 0.070 | −1.002 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A- T- P- | 0.030 | 0.076 | 0.399 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A- T- P- | −0.249 | 0.072 | −3.483 | 0.012 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T- P- | −0.332 | 0.091 | −3.665 | 0.006 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T- P- | −0.106 | 0.077 | −1.387 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T- P- | A- T- P+ | 0.092 | 0.072 | 1.276 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T- P+ | A- T- P+ | 0.000 | 0.070 | −0.005 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A- T- P+ | 0.100 | 0.075 | 1.326 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A- T- P+ | −0.180 | 0.071 | −2.525 | 0.182 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T- P+ | −0.262 | 0.090 | −2.906 | 0.070 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T- P+ | −0.037 | 0.076 | −0.480 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T- P+ | A+ T- P- | −0.093 | 0.073 | −1.265 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A+ T- P- | 0.008 | 0.079 | 0.096 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A+ T- P- | −0.272 | 0.075 | −3.632 | 0.007 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A+ T- P- | −0.355 | 0.093 | −3.804 | 0.004 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A+ T- P- | −0.129 | 0.080 | −1.617 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A+ T- P+ | 0.100 | 0.076 | 1.315 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A+ T- P+ | −0.179 | 0.072 | −2.488 | 0.198 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A+ T- P+ | −0.262 | 0.091 | −2.879 | 0.075 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A+ T- P+ | −0.036 | 0.077 | −0.470 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A- T+ P- | −0.280 | 0.078 | −3.594 | 0.007 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T+ P- | −0.362 | 0.096 | −3.790 | 0.004 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T+ P- | −0.137 | 0.083 | −1.655 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T+ P+ | −0.083 | 0.092 | −0.895 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T+ P+ | 0.143 | 0.079 | 1.818 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A+ T+ P- | 0.226 | 0.096 | 2.343 | ≥0.25 | Treatment |
| Richness | ||||||||
| Treatment | Alpha.Score | A- T- P+ | A- T- P- | 0.025 | 0.045 | 0.552 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T- P- | A- T- P- | −0.064 | 0.047 | −1.349 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T- P+ | A- T- P- | −0.108 | 0.046 | −2.363 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A- T- P- | −0.089 | 0.049 | −1.804 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A- T- P- | 0.022 | 0.047 | 0.462 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T- P- | 0.067 | 0.059 | 1.141 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T- P- | −0.028 | 0.050 | −0.560 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T- P- | A- T- P+ | −0.089 | 0.047 | −1.884 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T- P+ | A- T- P+ | −0.132 | 0.045 | −2.927 | 0.066 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A- T- P+ | −0.114 | 0.049 | −2.320 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A- T- P+ | −0.003 | 0.046 | −0.070 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T- P+ | 0.043 | 0.059 | 0.724 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T- P+ | −0.053 | 0.050 | −1.062 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T- P+ | A+ T- P- | −0.044 | 0.048 | −0.913 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A+ T- P- | −0.025 | 0.051 | −0.488 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A+ T- P- | 0.086 | 0.049 | 1.755 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A+ T- P- | 0.131 | 0.061 | 2.164 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A+ T- P- | 0.036 | 0.052 | 0.694 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A+ T- P+ | 0.019 | 0.050 | 0.373 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A+ T- P+ | 0.129 | 0.047 | 2.752 | 0.106 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A+ T- P+ | 0.175 | 0.059 | 2.954 | 0.061 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A+ T- P+ | 0.080 | 0.050 | 1.585 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A- T+ P- | 0.111 | 0.051 | 2.184 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T+ P- | 0.156 | 0.062 | 2.514 | 0.187 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T+ P- | 0.061 | 0.054 | 1.137 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T+ P+ | 0.046 | 0.060 | 0.762 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T+ P+ | −0.050 | 0.051 | −0.967 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A+ T+ P- | −0.095 | 0.063 | −1.520 | ≥0.25 | Treatment |
| Significant Results:Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE | ||||||||
| Tukey(Displacement ~ Treatment); All treatments at 60 DPE | ||||||||
| term | .y. | group1 | group2 | estimate | std.error | statistic | adj.p.value | Variable |
|---|---|---|---|---|---|---|---|---|
| Simpson | ||||||||
| Treatment | Alpha.Score | A- T+ P+ | A- T- P- | −0.249 | 0.072 | −3.483 | 0.012 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T- P- | −0.332 | 0.091 | −3.665 | 0.006 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A+ T- P- | −0.272 | 0.075 | −3.632 | 0.007 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A+ T- P- | −0.355 | 0.093 | −3.804 | 0.004 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A- T+ P- | −0.280 | 0.078 | −3.594 | 0.007 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T+ P- | −0.362 | 0.096 | −3.790 | 0.004 | Treatment |
# Set the subsection for this analysis
tmp.resSubSection <- "All"
# GLM Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]] <-
test.beta.all_metrics$volatility %>%
# Clean up cell value names by removing any strings including and after "__"
cutCellNames(col = "Beta.Metric", sep = "__") %>%
# Filter out raw data, leave only normalized data
dplyr::filter(Data.Type != "Raw") %>%
# Run GLM models
run_glm_models_beta(formula_str = "Displacement ~ Treatment",
beta_metric_col = "Beta.Metric",
beta_score_col = "Displacement",
family_str = "gaussian")
# ANOVA Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["ANOVA"]] <-
run_glm_anova_beta(beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]])
# Tukey Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["Tukey"]] <-
test.beta.all_metrics$volatility %>%
# Clean up cell value names by removing any strings including and after "__"
cutCellNames(col = "Beta.Metric", sep = "__") %>%
# Filter out raw data, leave only normalized data
dplyr::filter(Data.Type != "Raw") %>%
# Run Tukey test on GLM models
run_tukey_glm_beta(., "Displacement", "Beta.Metric", c("Treatment"),
family_str = "gaussian",
group_by_var = NULL)
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
| GLM Results - Beta | |||||
| glm(Displacement ~ Treatment); All treatments | |||||
| term | estimate | std.error | statistic | p.value | p.adj.sig |
|---|---|---|---|---|---|
| Bray-Curtis_norm | |||||
| (Intercept) | −0.061 | 0.027 | −2.266 | 0.024 | * |
| A- T- P+ | 0.142 | 0.038 | 3.716 | <0.001 | *** |
| A+ T- P- | 0.016 | 0.039 | 0.412 | ≥0.25 | ns |
| A+ T- P+ | 0.002 | 0.039 | 0.040 | ≥0.25 | ns |
| A- T+ P- | 0.081 | 0.040 | 2.048 | 0.041 | * |
| A- T+ P+ | −0.163 | 0.039 | −4.216 | <0.001 | **** |
| A+ T+ P- | −0.020 | 0.041 | −0.473 | ≥0.25 | ns |
| A+ T+ P+ | −0.004 | 0.040 | −0.104 | ≥0.25 | ns |
| Canberra_norm | |||||
| (Intercept) | −0.176 | 0.021 | −8.202 | <0.001 | **** |
| A- T- P+ | 0.070 | 0.030 | 2.300 | 0.022 | * |
| A+ T- P- | 0.040 | 0.031 | 1.307 | 0.192 | ns |
| A+ T- P+ | 0.083 | 0.031 | 2.726 | 0.007 | ** |
| A- T+ P- | 0.179 | 0.031 | 5.696 | <0.001 | **** |
| A- T+ P+ | 0.054 | 0.031 | 1.742 | 0.082 | ns |
| A+ T+ P- | 0.155 | 0.033 | 4.698 | <0.001 | **** |
| A+ T+ P+ | 0.096 | 0.032 | 3.044 | 0.002 | ** |
| ANOVA of GLM - Beta | ||||
| ANOVA(GLM(Displacement ~ Treatment), type = 2); All treatments | ||||
| term | statistic | df | p.value | sig |
|---|---|---|---|---|
| Bray-Curtis_norm | ||||
| Treatment | 70.305 | 7.000 | <0.001 | **** |
| Canberra_norm | ||||
| Treatment | 46.734 | 7.000 | <0.001 | **** |
| Pairwise Tukey's HSD, p.adj: Dunnett - Beta | ||||||||
| Tukey(Displacement ~ Treatment); All treatments | ||||||||
| term | .y. | group1 | group2 | estimate | std.error | statistic | adj.p.value | Variable |
|---|---|---|---|---|---|---|---|---|
| Bray-Curtis_norm | ||||||||
| Treatment | Beta.Score | A- T- P+ | A- T- P- | 0.142 | 0.038 | 3.716 | 0.005 | Treatment |
| Treatment | Beta.Score | A+ T- P- | A- T- P- | 0.016 | 0.039 | 0.412 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T- P+ | A- T- P- | 0.002 | 0.039 | 0.040 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P- | A- T- P- | 0.081 | 0.040 | 2.048 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A- T- P- | −0.163 | 0.039 | −4.216 | <0.001 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A- T- P- | −0.020 | 0.041 | −0.473 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A- T- P- | −0.004 | 0.040 | −0.104 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T- P- | A- T- P+ | −0.126 | 0.039 | −3.232 | 0.027 | Treatment |
| Treatment | Beta.Score | A+ T- P+ | A- T- P+ | −0.141 | 0.039 | −3.652 | 0.006 | Treatment |
| Treatment | Beta.Score | A- T+ P- | A- T- P+ | −0.061 | 0.040 | −1.543 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A- T- P+ | −0.306 | 0.039 | −7.884 | <0.001 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A- T- P+ | −0.162 | 0.041 | −3.900 | 0.002 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A- T- P+ | −0.146 | 0.040 | −3.682 | 0.006 | Treatment |
| Treatment | Beta.Score | A+ T- P+ | A+ T- P- | −0.015 | 0.039 | −0.370 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P- | A+ T- P- | 0.065 | 0.040 | 1.612 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A+ T- P- | −0.180 | 0.040 | −4.543 | <0.001 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A+ T- P- | −0.036 | 0.042 | −0.846 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A+ T- P- | −0.020 | 0.040 | −0.499 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P- | A+ T- P+ | 0.080 | 0.040 | 1.997 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A+ T- P+ | −0.165 | 0.039 | −4.230 | <0.001 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A+ T- P+ | −0.021 | 0.042 | −0.507 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A+ T- P+ | −0.006 | 0.040 | −0.142 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A- T+ P- | −0.245 | 0.040 | −6.102 | <0.001 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A- T+ P- | −0.101 | 0.043 | −2.358 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A- T+ P- | −0.085 | 0.041 | −2.077 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A- T+ P+ | 0.144 | 0.042 | 3.429 | 0.014 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A- T+ P+ | 0.159 | 0.040 | 3.961 | 0.002 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A+ T+ P- | 0.015 | 0.043 | 0.361 | ≥0.25 | Treatment |
| Canberra_norm | ||||||||
| Treatment | Beta.Score | A- T- P+ | A- T- P- | 0.070 | 0.030 | 2.300 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T- P- | A- T- P- | 0.040 | 0.031 | 1.307 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T- P+ | A- T- P- | 0.083 | 0.031 | 2.726 | 0.115 | Treatment |
| Treatment | Beta.Score | A- T+ P- | A- T- P- | 0.179 | 0.031 | 5.696 | <0.001 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A- T- P- | 0.054 | 0.031 | 1.742 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A- T- P- | 0.155 | 0.033 | 4.698 | <0.001 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A- T- P- | 0.096 | 0.032 | 3.044 | 0.048 | Treatment |
| Treatment | Beta.Score | A+ T- P- | A- T- P+ | −0.029 | 0.031 | −0.948 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T- P+ | A- T- P+ | 0.013 | 0.031 | 0.441 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P- | A- T- P+ | 0.109 | 0.031 | 3.473 | 0.012 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A- T- P+ | −0.016 | 0.031 | −0.528 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A- T- P+ | 0.085 | 0.033 | 2.576 | 0.164 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A- T- P+ | 0.026 | 0.032 | 0.830 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T- P+ | A+ T- P- | 0.043 | 0.031 | 1.374 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P- | A+ T- P- | 0.138 | 0.032 | 4.327 | <0.001 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A+ T- P- | 0.013 | 0.031 | 0.418 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A+ T- P- | 0.114 | 0.033 | 3.410 | 0.015 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A+ T- P- | 0.055 | 0.032 | 1.729 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P- | A+ T- P+ | 0.096 | 0.032 | 3.027 | 0.050 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A+ T- P+ | −0.030 | 0.031 | −0.960 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A+ T- P+ | 0.071 | 0.033 | 2.156 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A+ T- P+ | 0.013 | 0.032 | 0.400 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A- T+ P- | −0.125 | 0.032 | −3.943 | 0.002 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A- T+ P- | −0.024 | 0.034 | −0.718 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A- T+ P- | −0.083 | 0.033 | −2.549 | 0.175 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A- T+ P+ | 0.101 | 0.033 | 3.036 | 0.049 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A- T+ P+ | 0.042 | 0.032 | 1.329 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A+ T+ P- | −0.059 | 0.034 | −1.725 | ≥0.25 | Treatment |
# Set the subsection for this analysis
tmp.resSubSection <- "All_60DPE"
# GLM Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]] <-
test.beta.all_metrics$displacement %>%
dplyr::filter(Time == 60) %>%
# Clean up cell value names by removing any strings including and after "__"
cutCellNames(col = "Beta.Metric", sep = "__") %>%
# Filter out raw data, leave only normalized data
dplyr::filter(Data.Type != "Raw") %>%
# Run GLM models
run_glm_models_beta(formula_str = "Displacement ~ Treatment",
beta_metric_col = "Beta.Metric",
beta_score_col = "Displacement",
family_str = "gaussian")
# ANOVA Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["ANOVA"]] <-
run_glm_anova_beta(beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]])
# Tukey Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["Tukey"]] <-
test.beta.all_metrics$displacement %>%
dplyr::filter(Time == 60) %>%
# Clean up cell value names by removing any strings including and after "__"
cutCellNames(col = "Beta.Metric", sep = "__") %>%
# Filter out raw data, leave only normalized data
dplyr::filter(Data.Type != "Raw") %>%
# Run Tukey test on GLM models
run_tukey_glm_beta(., "Displacement", "Beta.Metric", c("Treatment"),
family_str = "gaussian",
group_by_var = NULL)
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
| GLM Results - 60 DPE | |||||
| glm(Displacement ~ Treatment); All treatments at 60 DPE | |||||
| term | estimate | std.error | statistic | p.value | p.adj.sig |
|---|---|---|---|---|---|
| Bray-Curtis_norm | |||||
| (Intercept) | −0.060 | 0.041 | −1.474 | 0.142 | ns |
| A- T- P+ | 0.284 | 0.057 | 4.945 | <0.001 | **** |
| A+ T- P- | 0.007 | 0.061 | 0.115 | ≥0.25 | ns |
| A+ T- P+ | −0.021 | 0.058 | −0.360 | ≥0.25 | ns |
| A- T+ P- | 0.043 | 0.063 | 0.675 | ≥0.25 | ns |
| A- T+ P+ | −0.206 | 0.060 | −3.458 | <0.001 | *** |
| A+ T+ P- | −0.119 | 0.075 | −1.584 | 0.115 | ns |
| A+ T+ P+ | 0.063 | 0.064 | 0.993 | ≥0.25 | ns |
| Canberra_norm | |||||
| (Intercept) | −0.124 | 0.034 | −3.647 | <0.001 | *** |
| A- T- P+ | 0.134 | 0.048 | 2.796 | 0.006 | ** |
| A+ T- P- | −0.089 | 0.051 | −1.750 | 0.081 | ns |
| A+ T- P+ | 0.037 | 0.049 | 0.765 | ≥0.25 | ns |
| A- T+ P- | 0.099 | 0.053 | 1.871 | 0.063 | ns |
| A- T+ P+ | 0.065 | 0.050 | 1.306 | 0.193 | ns |
| A+ T+ P- | 0.123 | 0.063 | 1.962 | 0.051 | ns |
| A+ T+ P+ | 0.071 | 0.053 | 1.329 | 0.185 | ns |
| ANOVA of GLM - 60 DPE | ||||
| ANOVA(GLM(Displacement ~ Treatment), type = 2); All treatments at 60 DPE | ||||
| term | statistic | df | p.value | sig |
|---|---|---|---|---|
| Bray-Curtis_norm | ||||
| Treatment | 77.159 | 7.000 | <0.001 | **** |
| Canberra_norm | ||||
| Treatment | 26.083 | 7.000 | <0.001 | *** |
| Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE | ||||||||
| Tukey(Displacement ~ Treatment); All treatments at 60 DPE | ||||||||
| term | .y. | group1 | group2 | estimate | std.error | statistic | adj.p.value | Variable |
|---|---|---|---|---|---|---|---|---|
| Bray-Curtis_norm | ||||||||
| Treatment | Beta.Score | A- T- P+ | A- T- P- | 0.284 | 0.057 | 4.945 | <0.001 | Treatment |
| Treatment | Beta.Score | A+ T- P- | A- T- P- | 0.007 | 0.061 | 0.115 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T- P+ | A- T- P- | −0.021 | 0.058 | −0.360 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P- | A- T- P- | 0.043 | 0.063 | 0.675 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A- T- P- | −0.206 | 0.060 | −3.458 | 0.013 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A- T- P- | −0.119 | 0.075 | −1.584 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A- T- P- | 0.063 | 0.064 | 0.993 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T- P- | A- T- P+ | −0.277 | 0.060 | −4.595 | <0.001 | Treatment |
| Treatment | Beta.Score | A+ T- P+ | A- T- P+ | −0.305 | 0.058 | −5.272 | <0.001 | Treatment |
| Treatment | Beta.Score | A- T+ P- | A- T- P+ | −0.241 | 0.063 | −3.845 | 0.003 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A- T- P+ | −0.490 | 0.059 | −8.275 | <0.001 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A- T- P+ | −0.403 | 0.075 | −5.372 | <0.001 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A- T- P+ | −0.220 | 0.063 | −3.473 | 0.012 | Treatment |
| Treatment | Beta.Score | A+ T- P+ | A+ T- P- | −0.028 | 0.061 | −0.458 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P- | A+ T- P- | 0.036 | 0.066 | 0.542 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A+ T- P- | −0.213 | 0.062 | −3.418 | 0.014 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A+ T- P- | −0.126 | 0.078 | −1.629 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A+ T- P- | 0.056 | 0.066 | 0.849 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P- | A+ T- P+ | 0.064 | 0.063 | 1.001 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A+ T- P+ | −0.185 | 0.060 | −3.086 | 0.042 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A+ T- P+ | −0.098 | 0.076 | −1.300 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A+ T- P+ | 0.084 | 0.064 | 1.313 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A- T+ P- | −0.249 | 0.065 | −3.840 | 0.003 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A- T+ P- | −0.162 | 0.079 | −2.037 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A- T+ P- | 0.021 | 0.069 | 0.303 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A- T+ P+ | 0.087 | 0.077 | 1.130 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A- T+ P+ | 0.269 | 0.065 | 4.116 | <0.001 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A+ T+ P- | 0.183 | 0.080 | 2.282 | ≥0.25 | Treatment |
| Canberra_norm | ||||||||
| Treatment | Beta.Score | A- T- P+ | A- T- P- | 0.134 | 0.048 | 2.796 | 0.094 | Treatment |
| Treatment | Beta.Score | A+ T- P- | A- T- P- | −0.089 | 0.051 | −1.750 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T- P+ | A- T- P- | 0.037 | 0.049 | 0.765 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P- | A- T- P- | 0.099 | 0.053 | 1.871 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A- T- P- | 0.065 | 0.050 | 1.306 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A- T- P- | 0.123 | 0.063 | 1.962 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A- T- P- | 0.071 | 0.053 | 1.329 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T- P- | A- T- P+ | −0.223 | 0.050 | −4.425 | <0.001 | Treatment |
| Treatment | Beta.Score | A+ T- P+ | A- T- P+ | −0.097 | 0.048 | −2.006 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P- | A- T- P+ | −0.035 | 0.052 | −0.676 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A- T- P+ | −0.069 | 0.049 | −1.397 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A- T- P+ | −0.011 | 0.063 | −0.168 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A- T- P+ | −0.063 | 0.053 | −1.192 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T- P+ | A+ T- P- | 0.126 | 0.051 | 2.468 | 0.206 | Treatment |
| Treatment | Beta.Score | A- T+ P- | A+ T- P- | 0.187 | 0.055 | 3.412 | 0.014 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A+ T- P- | 0.153 | 0.052 | 2.951 | 0.062 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A+ T- P- | 0.212 | 0.065 | 3.275 | 0.023 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A+ T- P- | 0.159 | 0.055 | 2.875 | 0.076 | Treatment |
| Treatment | Beta.Score | A- T+ P- | A+ T- P+ | 0.061 | 0.053 | 1.159 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A+ T- P+ | 0.028 | 0.050 | 0.555 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A+ T- P+ | 0.086 | 0.063 | 1.366 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A+ T- P+ | 0.034 | 0.054 | 0.628 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A- T+ P- | −0.034 | 0.054 | −0.623 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A- T+ P- | 0.025 | 0.066 | 0.375 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A- T+ P- | −0.028 | 0.057 | −0.485 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A- T+ P+ | 0.058 | 0.064 | 0.913 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A- T+ P+ | 0.006 | 0.055 | 0.107 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A+ T+ P- | −0.053 | 0.067 | −0.788 | ≥0.25 | Treatment |
| Significant Results:Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE | ||||||||
| Tukey(Displacement ~ Treatment); All treatments at 60 DPE | ||||||||
| term | .y. | group1 | group2 | estimate | std.error | statistic | adj.p.value | Variable |
|---|---|---|---|---|---|---|---|---|
| Bray-Curtis_norm | ||||||||
| Treatment | Beta.Score | A- T- P+ | A- T- P- | 0.284 | 0.057 | 4.945 | <0.001 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A- T- P- | −0.206 | 0.060 | −3.458 | 0.013 | Treatment |
| Treatment | Beta.Score | A+ T- P- | A- T- P+ | −0.277 | 0.060 | −4.595 | <0.001 | Treatment |
| Treatment | Beta.Score | A+ T- P+ | A- T- P+ | −0.305 | 0.058 | −5.272 | <0.001 | Treatment |
| Treatment | Beta.Score | A- T+ P- | A- T- P+ | −0.241 | 0.063 | −3.845 | 0.003 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A- T- P+ | −0.490 | 0.059 | −8.275 | <0.001 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A- T- P+ | −0.403 | 0.075 | −5.372 | <0.001 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A- T- P+ | −0.220 | 0.063 | −3.473 | 0.012 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A+ T- P- | −0.213 | 0.062 | −3.418 | 0.014 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A+ T- P+ | −0.185 | 0.060 | −3.086 | 0.042 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A- T+ P- | −0.249 | 0.065 | −3.840 | 0.003 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A- T+ P+ | 0.269 | 0.065 | 4.116 | <0.001 | Treatment |
| Canberra_norm | ||||||||
| Treatment | Beta.Score | A+ T- P- | A- T- P+ | −0.223 | 0.050 | −4.425 | <0.001 | Treatment |
| Treatment | Beta.Score | A- T+ P- | A+ T- P- | 0.187 | 0.055 | 3.412 | 0.014 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A+ T- P- | 0.212 | 0.065 | 3.275 | 0.023 | Treatment |
# Set the subsection for this analysis
tmp.resSubSection <- "All_60DPE"
#### Dispersion --------------------------------------------------------------
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]] <-
run_BetaDispersion(dist.matrix = beta.dist.mat[[tmp.resSubSection]],
var = c("treatment_group"))
##### ANOVA ------------------------------------------------------------------
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA"]] <-
get_HoD_anova(betaDisper = beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]],
var = c("treatment_group"))
### Table
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA.Table"]] <-
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA"]] %>%
# Create the table
dplyr::group_by(Beta.Metric) %>%
set_GT(var = "p.value", group.by = "Beta.Metric") %>%
# Title/caption
gt::tab_header(
title = "ANOVA: Homogeneity of Dispersion",
subtitle = "ANOVA(Beta Disperson ~ Treatment); All treatments at 60 DPE"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "lightgreen"),
locations = gt::cells_body(
columns = c("p.value", "sig"),
rows = p.value < 0.05
)
)
##### Tukey ------------------------------------------------------------------
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey"]] <-
get_HoD_tukey(betaDisper = beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]],
var = c("treatment_group"))
### Table
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey.Table"]] <-
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey"]] %>%
dplyr::mutate(
group1 = case_when(
group1 == "Control" ~ "A- T- P-",
group1 == "Parasite" ~ "A- T- P+",
group1 == "Antibiotics" ~ "A+ T- P-",
group1 == "Antibiotics_Parasite" ~ "A+ T- P+",
group1 == "Temperature" ~ "A- T+ P-",
group1 == "Temperature_Parasite" ~ "A- T+ P+",
group1 == "Antibiotics_Temperature" ~ "A+ T+ P-",
group1 == "Antibiotics_Temperature_Parasite" ~ "A+ T+ P+",
),
group2 = case_when(
group2 == "Control" ~ "A- T- P-",
group2 == "Parasite" ~ "A- T- P+",
group2 == "Antibiotics" ~ "A+ T- P-",
group2 == "Antibiotics_Parasite" ~ "A+ T- P+",
group2 == "Temperature" ~ "A- T+ P-",
group2 == "Temperature_Parasite" ~ "A- T+ P+",
group2 == "Antibiotics_Temperature" ~ "A+ T+ P-",
group2 == "Antibiotics_Temperature_Parasite" ~ "A+ T+ P+",
)) %>%
# Create the table
dplyr::group_by(Beta.Metric) %>%
set_GT(var = "adj.p.value", group.by = "Beta.Metric") %>%
# Title/caption
gt::tab_header(
title = "Tukey: Homogeneity of Dispersion",
subtitle = "Tukey(Beta Disperson ~ Treatment); All treatments at 60 DPE"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "lightgreen"),
locations = gt::cells_body(
columns = c("adj.p.value", "sig"),
rows = adj.p.value < 0.05
)
)
| ANOVA: Homogeneity of Dispersion | ||||||
| ANOVA(Beta Disperson ~ Treatment); All treatments at 60 DPE | ||||||
| term | df | sumsq | meansq | statistic | p.value | sig |
|---|---|---|---|---|---|---|
| bray | ||||||
| treatment_group | 7.000 | 0.645 | 0.092 | 3.523 | 0.001 | ** |
| Residual | 228.000 | 5.965 | 0.026 | NA | NA | NA |
| canberra | ||||||
| treatment_group | 7.000 | 0.132 | 0.019 | 5.021 | <0.001 | **** |
| Residual | 228.000 | 0.857 | 0.004 | NA | NA | NA |
| Tukey: Homogeneity of Dispersion | ||||||||
| Tukey(Beta Disperson ~ Treatment); All treatments at 60 DPE | ||||||||
| .y. | term | group1 | group2 | estimate | conf.low | conf.high | adj.p.value | sig |
|---|---|---|---|---|---|---|---|---|
| bray | ||||||||
| Distance | treatment_group | A- T- P+ | A- T- P- | −0.025 | −0.141 | 0.091 | ≥0.25 | ns |
| Distance | treatment_group | A+ T- P- | A- T- P- | −0.036 | −0.158 | 0.086 | ≥0.25 | ns |
| Distance | treatment_group | A+ T- P+ | A- T- P- | −0.065 | −0.183 | 0.052 | ≥0.25 | ns |
| Distance | treatment_group | A- T+ P- | A- T- P- | −0.064 | −0.191 | 0.064 | ≥0.25 | ns |
| Distance | treatment_group | A- T+ P+ | A- T- P- | −0.106 | −0.227 | 0.014 | 0.126 | ns |
| Distance | treatment_group | A+ T+ P- | A- T- P- | −0.162 | −0.314 | −0.010 | 0.028 | * |
| Distance | treatment_group | A+ T+ P+ | A- T- P- | −0.157 | −0.286 | −0.028 | 0.006 | ** |
| Distance | treatment_group | A+ T- P- | A- T- P+ | −0.011 | −0.132 | 0.111 | ≥0.25 | ns |
| Distance | treatment_group | A+ T- P+ | A- T- P+ | −0.040 | −0.157 | 0.076 | ≥0.25 | ns |
| Distance | treatment_group | A- T+ P- | A- T- P+ | −0.039 | −0.165 | 0.088 | ≥0.25 | ns |
| Distance | treatment_group | A- T+ P+ | A- T- P+ | −0.081 | −0.201 | 0.038 | ≥0.25 | ns |
| Distance | treatment_group | A+ T+ P- | A- T- P+ | −0.137 | −0.289 | 0.014 | 0.108 | ns |
| Distance | treatment_group | A+ T+ P+ | A- T- P+ | −0.132 | −0.260 | −0.004 | 0.038 | * |
| Distance | treatment_group | A+ T- P+ | A+ T- P- | −0.029 | −0.152 | 0.094 | ≥0.25 | ns |
| Distance | treatment_group | A- T+ P- | A+ T- P- | −0.028 | −0.160 | 0.105 | ≥0.25 | ns |
| Distance | treatment_group | A- T+ P+ | A+ T- P- | −0.070 | −0.196 | 0.055 | ≥0.25 | ns |
| Distance | treatment_group | A+ T+ P- | A+ T- P- | −0.126 | −0.283 | 0.030 | 0.215 | ns |
| Distance | treatment_group | A+ T+ P+ | A+ T- P- | −0.121 | −0.255 | 0.013 | 0.108 | ns |
| Distance | treatment_group | A- T+ P- | A+ T- P+ | 0.001 | −0.127 | 0.130 | ≥0.25 | ns |
| Distance | treatment_group | A- T+ P+ | A+ T- P+ | −0.041 | −0.162 | 0.080 | ≥0.25 | ns |
| Distance | treatment_group | A+ T+ P- | A+ T- P+ | −0.097 | −0.250 | 0.056 | ≥0.25 | ns |
| Distance | treatment_group | A+ T+ P+ | A+ T- P+ | −0.092 | −0.222 | 0.038 | ≥0.25 | ns |
| Distance | treatment_group | A- T+ P+ | A- T+ P- | −0.042 | −0.173 | 0.088 | ≥0.25 | ns |
| Distance | treatment_group | A+ T+ P- | A- T+ P- | −0.098 | −0.259 | 0.062 | ≥0.25 | ns |
| Distance | treatment_group | A+ T+ P+ | A- T+ P- | −0.093 | −0.232 | 0.045 | ≥0.25 | ns |
| Distance | treatment_group | A+ T+ P- | A- T+ P+ | −0.056 | −0.211 | 0.099 | ≥0.25 | ns |
| Distance | treatment_group | A+ T+ P+ | A- T+ P+ | −0.051 | −0.183 | 0.081 | ≥0.25 | ns |
| Distance | treatment_group | A+ T+ P+ | A+ T+ P- | 0.005 | −0.157 | 0.166 | ≥0.25 | ns |
| canberra | ||||||||
| Distance | treatment_group | A- T- P+ | A- T- P- | 0.008 | −0.036 | 0.052 | ≥0.25 | ns |
| Distance | treatment_group | A+ T- P- | A- T- P- | −0.056 | −0.103 | −0.010 | 0.006 | ** |
| Distance | treatment_group | A+ T- P+ | A- T- P- | 0.007 | −0.037 | 0.052 | ≥0.25 | ns |
| Distance | treatment_group | A- T+ P- | A- T- P- | 0.010 | −0.039 | 0.058 | ≥0.25 | ns |
| Distance | treatment_group | A- T+ P+ | A- T- P- | −0.007 | −0.052 | 0.039 | ≥0.25 | ns |
| Distance | treatment_group | A+ T+ P- | A- T- P- | −0.023 | −0.080 | 0.035 | ≥0.25 | ns |
| Distance | treatment_group | A+ T+ P+ | A- T- P- | −0.044 | −0.093 | 0.005 | 0.116 | ns |
| Distance | treatment_group | A+ T- P- | A- T- P+ | −0.065 | −0.111 | −0.019 | <0.001 | *** |
| Distance | treatment_group | A+ T- P+ | A- T- P+ | −0.001 | −0.045 | 0.043 | ≥0.25 | ns |
| Distance | treatment_group | A- T+ P- | A- T- P+ | 0.001 | −0.047 | 0.049 | ≥0.25 | ns |
| Distance | treatment_group | A- T+ P+ | A- T- P+ | −0.015 | −0.060 | 0.030 | ≥0.25 | ns |
| Distance | treatment_group | A+ T+ P- | A- T- P+ | −0.031 | −0.088 | 0.026 | ≥0.25 | ns |
| Distance | treatment_group | A+ T+ P+ | A- T- P+ | −0.052 | −0.101 | −0.004 | 0.026 | * |
| Distance | treatment_group | A+ T- P+ | A+ T- P- | 0.064 | 0.017 | 0.110 | 0.001 | ** |
| Distance | treatment_group | A- T+ P- | A+ T- P- | 0.066 | 0.016 | 0.116 | 0.002 | ** |
| Distance | treatment_group | A- T+ P+ | A+ T- P- | 0.050 | 0.002 | 0.097 | 0.034 | * |
| Distance | treatment_group | A+ T+ P- | A+ T- P- | 0.034 | −0.026 | 0.093 | ≥0.25 | ns |
| Distance | treatment_group | A+ T+ P+ | A+ T- P- | 0.013 | −0.038 | 0.063 | ≥0.25 | ns |
| Distance | treatment_group | A- T+ P- | A+ T- P+ | 0.002 | −0.046 | 0.051 | ≥0.25 | ns |
| Distance | treatment_group | A- T+ P+ | A+ T- P+ | −0.014 | −0.060 | 0.032 | ≥0.25 | ns |
| Distance | treatment_group | A+ T+ P- | A+ T- P+ | −0.030 | −0.088 | 0.028 | ≥0.25 | ns |
| Distance | treatment_group | A+ T+ P+ | A+ T- P+ | −0.051 | −0.100 | −0.002 | 0.035 | * |
| Distance | treatment_group | A- T+ P+ | A- T+ P- | −0.016 | −0.066 | 0.033 | ≥0.25 | ns |
| Distance | treatment_group | A+ T+ P- | A- T+ P- | −0.032 | −0.093 | 0.029 | ≥0.25 | ns |
| Distance | treatment_group | A+ T+ P+ | A- T+ P- | −0.053 | −0.106 | −0.001 | 0.044 | * |
| Distance | treatment_group | A+ T+ P- | A- T+ P+ | −0.016 | −0.075 | 0.043 | ≥0.25 | ns |
| Distance | treatment_group | A+ T+ P+ | A- T+ P+ | −0.037 | −0.087 | 0.013 | ≥0.25 | ns |
| Distance | treatment_group | A+ T+ P+ | A+ T+ P- | −0.021 | −0.082 | 0.040 | ≥0.25 | ns |
# Set the subsection for this analysis
tmp.resSubSection <- "All_60DPE"
#### Capscale ----------------------------------------------------------------
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.mod"]] <-
run_capscale(ps.cleaned_alpha %>% # is the phyloseq object we are looping over
# microViz::ps_filter(Treatment %in% c(
# "A- T- P-",
# "A- T- P+",
# "A+ T- P-",
# "A+ T- P+",
# "A- T+ P-",
# "A- T+ P+",
# "A+ T+ P-",
# "A+ T+ P+"
# )) %>%
microViz::ps_filter(Time == 60),
dist.matrix = beta.dist.mat[[tmp.resSubSection]],
formula_str = "dist ~ Treatment")
##### ADONIS ------------------------------------------------------------------
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS"]] <-
run_cap_adonis(ps.cleaned_alpha %>% # is the phyloseq object we are looping over
# microViz::ps_filter(Treatment %in% c(
# "A- T- P-",
# "A- T- P+",
# "A+ T- P-",
# "A+ T- P+",
# "A- T+ P-",
# "A- T+ P+",
# "A+ T+ P-",
# "A+ T+ P+"
# )) %>%
microViz::ps_filter(Time == 60),
dist.matrix = beta.dist.mat[[tmp.resSubSection]],
formula_str = "dist ~ Treatment",
by.method = "terms")
## Warning in tidy.anova(.): The following column names in ANOVA output were not
## recognized or transformed: SumOfSqs, R2
## Warning in tidy.anova(.): The following column names in ANOVA output were not
## recognized or transformed: SumOfSqs, R2
### Table
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS.Table"]] <-
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS"]] %>%
set_GT(var = "p.value", group.by = "Beta.Metric") %>%
# Title/caption
gt::tab_header(
title = "ADONIS2",
subtitle = "adonis2(Beta Distance ~ Treatment); All treatments at 60 DPE"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "lightgreen"),
locations = gt::cells_body(
columns = c("p.value", "sig"),
rows = p.value < 0.05
)
)
| ADONIS2 | ||||||
| adonis2(Beta Distance ~ Treatment); All treatments at 60 DPE | ||||||
| term | df | SumOfSqs | R2 | statistic | p.value | sig |
|---|---|---|---|---|---|---|
| bray | ||||||
| Treatment | 7.000 | 3.508 | 0.105 | 3.813 | 0.001 | *** |
| Residual | 228.000 | 29.963 | 0.895 | NA | NA | NA |
| Total | 235.000 | 33.471 | 1.000 | NA | NA | NA |
| canberra | ||||||
| Treatment | 7.000 | 8.648 | 0.103 | 3.733 | 0.001 | *** |
| Residual | 228.000 | 75.464 | 0.897 | NA | NA | NA |
| Total | 235.000 | 84.112 | 1.000 | NA | NA | NA |
# Set the subsection for this analysis
tmp.resSubSection <- "ParasiteExp"
# GLM Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]] <-
test.alpha.all_metrics$volatility %>%
# Filter
dplyr::filter(Treatment %in% c(
"A- T- P-",
"A- T- P+"#,
# "A+ T- P-",
# "A+ T- P+",
# "A- T+ P-",
# "A- T+ P+",
# "A+ T+ P-",
# "A+ T+ P+"
)) %>%
# Clean up cell value names by removing any strings including and after "__"
cutCellNames(col = "Alpha.Metric", sep = "__") %>%
# Filter out raw data, leave only normalized data
dplyr::filter(Data.Type != "Raw") %>%
# Run GLM models
run_glm_models(formula_str = "Displacement ~ Treatment",
alpha_metric_col = "Alpha.Metric",
alpha_score_col = "Displacement",
family_str = "gaussian")
# ANOVA Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["ANOVA"]] <-
run_glm_anova(alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]])
# Tukey Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["Tukey"]] <-
test.alpha.all_metrics$volatility %>%
# Filter
dplyr::filter(Treatment %in% c(
"A- T- P-",
"A- T- P+"#,
# "A+ T- P-",
# "A+ T- P+",
# "A- T+ P-",
# "A- T+ P+",
# "A+ T+ P-",
# "A+ T+ P+"
)) %>%
# Clean up cell value names by removing any strings including and after "__"
cutCellNames(col = "Alpha.Metric", sep = "__") %>%
# Filter out raw data, leave only normalized data
dplyr::filter(Data.Type != "Raw") %>%
# Run Tukey test on GLM models
run_tukey_glm(., "Displacement", "Alpha.Metric", c("Treatment"),
family_str = "gaussian",
group_by_var = NULL)
| GLM Results | |||||
| glm(Displacement ~ Treatment); Parasite exposure only | |||||
| term | estimate | std.error | statistic | p.value | p.adj.sig |
|---|---|---|---|---|---|
| Shannon | |||||
| (Intercept) | 0.078 | 0.021 | 3.755 | <0.001 | *** |
| A- T- P+ | −0.045 | 0.029 | −1.529 | 0.128 | ns |
| Simpson | |||||
| (Intercept) | 0.042 | 0.029 | 1.480 | 0.141 | ns |
| A- T- P+ | −0.049 | 0.040 | −1.199 | 0.232 | ns |
| Richness | |||||
| (Intercept) | 0.221 | 0.025 | 8.745 | <0.001 | **** |
| A- T- P+ | 0.009 | 0.036 | 0.251 | ≥0.25 | ns |
| ANOVA of GLM | ||||
| ANOVA(GLM(Displacement ~ Treatment), type = 2); Parasite exposure only | ||||
| term | statistic | df | p.value | sig |
|---|---|---|---|---|
| Shannon | ||||
| Treatment | 2.338 | 1.000 | 0.126 | ns |
| Simpson | ||||
| Treatment | 1.439 | 1.000 | 0.230 | ns |
| Richness | ||||
| Treatment | 0.063 | 1.000 | ≥0.25 | ns |
| Pairwise Tukey's HSD, p.adj: Dunnett | ||||||||
| Tukey(Displacement ~ Treatment); Parasite exposure only | ||||||||
| term | .y. | group1 | group2 | estimate | std.error | statistic | adj.p.value | Variable |
|---|---|---|---|---|---|---|---|---|
| Shannon | ||||||||
| Treatment | Alpha.Score | A- T- P+ | A- T- P- | −0.045 | 0.029 | −1.529 | 0.126 | Treatment |
| Simpson | ||||||||
| Treatment | Alpha.Score | A- T- P+ | A- T- P- | −0.049 | 0.040 | −1.199 | 0.230 | Treatment |
| Richness | ||||||||
| Treatment | Alpha.Score | A- T- P+ | A- T- P- | 0.009 | 0.036 | 0.251 | ≥0.25 | Treatment |
# Set the subsection for this analysis
tmp.resSubSection <- "ParasiteExp_60DPE"
# GLM Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]] <-
test.alpha.all_metrics$displacement %>%
dplyr::filter(Treatment %in% c(
"A- T- P-",
"A- T- P+"#,
# "A+ T- P-",
# "A+ T- P+",
# "A- T+ P-",
# "A- T+ P+",
# "A+ T+ P-",
# "A+ T+ P+"
)) %>%
dplyr::filter(Time == 60) %>%
# Clean up cell value names by removing any strings including and after "__"
cutCellNames(col = "Alpha.Metric", sep = "__") %>%
# Filter out raw data, leave only normalized data
dplyr::filter(Data.Type != "Raw") %>%
# Run GLM models
run_glm_models(formula_str = "Displacement ~ Treatment",
alpha_metric_col = "Alpha.Metric",
alpha_score_col = "Displacement",
family_str = "gaussian")
# ANOVA Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["ANOVA"]] <-
run_glm_anova(alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]])
# Tukey Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["Tukey"]] <-
test.alpha.all_metrics$displacement %>%
dplyr::filter(Treatment %in% c(
"A- T- P-",
"A- T- P+"#,
# "A+ T- P-",
# "A+ T- P+",
# "A- T+ P-",
# "A- T+ P+",
# "A+ T+ P-",
# "A+ T+ P+"
)) %>%
dplyr::filter(Time == 60) %>%
# Clean up cell value names by removing any strings including and after "__"
cutCellNames(col = "Alpha.Metric", sep = "__") %>%
# Filter out raw data, leave only normalized data
dplyr::filter(Data.Type != "Raw") %>%
# Run Tukey test on GLM models
run_tukey_glm(., "Displacement", "Alpha.Metric", c("Treatment"),
family_str = "gaussian",
group_by_var = NULL)
| GLM Results - 60 DPE | |||||
| glm(Displacement ~ Treatment); All treatments at 60 DPE | |||||
| term | estimate | std.error | statistic | p.value | p.adj.sig |
|---|---|---|---|---|---|
| Shannon | |||||
| (Intercept) | 0.091 | 0.041 | 2.207 | 0.031 | * |
| A- T- P+ | −0.082 | 0.058 | −1.428 | 0.158 | ns |
| Simpson | |||||
| (Intercept) | 0.008 | 0.053 | 0.150 | ≥0.25 | ns |
| A- T- P+ | −0.070 | 0.075 | −0.933 | ≥0.25 | ns |
| Richness | |||||
| (Intercept) | 0.338 | 0.034 | 9.916 | <0.001 | **** |
| A- T- P+ | 0.025 | 0.048 | 0.518 | ≥0.25 | ns |
| ANOVA of GLM - 60 DPE | ||||
| ANOVA(GLM(Displacement ~ Treatment), type = 2); Parasite exposure only at 60 DPE | ||||
| term | statistic | df | p.value | sig |
|---|---|---|---|---|
| Shannon | ||||
| Treatment | 2.039 | 1.000 | 0.153 | ns |
| Simpson | ||||
| Treatment | 0.871 | 1.000 | ≥0.25 | ns |
| Richness | ||||
| Treatment | 0.269 | 1.000 | ≥0.25 | ns |
| Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE | ||||||||
| Tukey(Displacement ~ Treatment); Parasite exposure only at 60 DPE | ||||||||
| term | .y. | group1 | group2 | estimate | std.error | statistic | adj.p.value | Variable |
|---|---|---|---|---|---|---|---|---|
| Shannon | ||||||||
| Treatment | Alpha.Score | A- T- P+ | A- T- P- | −0.082 | 0.058 | −1.428 | 0.153 | Treatment |
| Simpson | ||||||||
| Treatment | Alpha.Score | A- T- P+ | A- T- P- | −0.070 | 0.075 | −0.933 | ≥0.25 | Treatment |
| Richness | ||||||||
| Treatment | Alpha.Score | A- T- P+ | A- T- P- | 0.025 | 0.048 | 0.518 | ≥0.25 | Treatment |
| Significant Results:Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE | ||||||||
| Tukey(Displacement ~ Treatment); All treatments at 60 DPE | ||||||||
| term | .y. | group1 | group2 | estimate | std.error | statistic | adj.p.value | Variable |
|---|---|---|---|---|---|---|---|---|
# Set the subsection for this analysis
tmp.resSubSection <- "ParasiteExp"
# GLM Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]] <-
test.beta.all_metrics$volatility %>%
# Filter
dplyr::filter(Treatment %in% c(
"A- T- P-",
"A- T- P+"#,
# "A+ T- P-",
# "A+ T- P+",
# "A- T+ P-",
# "A- T+ P+",
# "A+ T+ P-",
# "A+ T+ P+"
)) %>%
# Clean up cell value names by removing any strings including and after "__"
cutCellNames(col = "Beta.Metric", sep = "__") %>%
# Filter out raw data, leave only normalized data
dplyr::filter(Data.Type != "Raw") %>%
# Run GLM models
run_glm_models_beta(formula_str = "Displacement ~ Treatment",
beta_metric_col = "Beta.Metric",
beta_score_col = "Displacement",
family_str = "gaussian")
# ANOVA Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["ANOVA"]] <-
run_glm_anova_beta(beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]])
# Tukey Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["Tukey"]] <-
test.beta.all_metrics$volatility %>%
# Filter
dplyr::filter(Treatment %in% c(
"A- T- P-",
"A- T- P+"#,
# "A+ T- P-",
# "A+ T- P+",
# "A- T+ P-",
# "A- T+ P+",
# "A+ T+ P-",
# "A+ T+ P+"
)) %>%
# Clean up cell value names by removing any strings including and after "__"
cutCellNames(col = "Beta.Metric", sep = "__") %>%
# Filter out raw data, leave only normalized data
dplyr::filter(Data.Type != "Raw") %>%
# Run Tukey test on GLM models
run_tukey_glm_beta(., "Displacement", "Beta.Metric", c("Treatment"),
family_str = "gaussian",
group_by_var = NULL)
| GLM Results - Beta | |||||
| glm(Displacement ~ Treatment); Parasite exposure only | |||||
| term | estimate | std.error | statistic | p.value | p.adj.sig |
|---|---|---|---|---|---|
| Bray-Curtis_norm | |||||
| (Intercept) | −0.061 | 0.026 | −2.369 | 0.019 | * |
| A- T- P+ | 0.142 | 0.037 | 3.884 | <0.001 | *** |
| Canberra_norm | |||||
| (Intercept) | −0.176 | 0.023 | −7.710 | <0.001 | **** |
| A- T- P+ | 0.070 | 0.032 | 2.162 | 0.032 | * |
| ANOVA of GLM - Beta | ||||
| ANOVA(GLM(Displacement ~ Treatment), type = 2); Parasite exposure only | ||||
| term | statistic | df | p.value | sig |
|---|---|---|---|---|
| Bray-Curtis_norm | ||||
| Treatment | 15.083 | 1.000 | <0.001 | *** |
| Canberra_norm | ||||
| Treatment | 4.673 | 1.000 | 0.031 | * |
| Pairwise Tukey's HSD, p.adj: Dunnett - Beta | ||||||||
| Tukey(Displacement ~ Treatment); Parasite exposure only | ||||||||
| term | .y. | group1 | group2 | estimate | std.error | statistic | adj.p.value | Variable |
|---|---|---|---|---|---|---|---|---|
| Bray-Curtis_norm | ||||||||
| Treatment | Beta.Score | A- T- P+ | A- T- P- | 0.142 | 0.037 | 3.884 | <0.001 | Treatment |
| Canberra_norm | ||||||||
| Treatment | Beta.Score | A- T- P+ | A- T- P- | 0.070 | 0.032 | 2.162 | 0.031 | Treatment |
# Set the subsection for this analysis
tmp.resSubSection <- "ParasiteExp_60DPE"
# GLM Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]] <-
test.beta.all_metrics$displacement %>%
dplyr::filter(Treatment %in% c(
"A- T- P-",
"A- T- P+"#,
# "A+ T- P-",
# "A+ T- P+",
# "A- T+ P-",
# "A- T+ P+",
# "A+ T+ P-",
# "A+ T+ P+"
)) %>%
dplyr::filter(Time == 60) %>%
# Clean up cell value names by removing any strings including and after "__"
cutCellNames(col = "Beta.Metric", sep = "__") %>%
# Filter out raw data, leave only normalized data
dplyr::filter(Data.Type != "Raw") %>%
# Run GLM models
run_glm_models_beta(formula_str = "Displacement ~ Treatment",
beta_metric_col = "Beta.Metric",
beta_score_col = "Displacement",
family_str = "gaussian")
# ANOVA Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["ANOVA"]] <-
run_glm_anova_beta(beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]])
# Tukey Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["Tukey"]] <-
test.beta.all_metrics$displacement %>%
dplyr::filter(Treatment %in% c(
"A- T- P-",
"A- T- P+"#,
# "A+ T- P-",
# "A+ T- P+",
# "A- T+ P-",
# "A- T+ P+",
# "A+ T+ P-",
# "A+ T+ P+"
)) %>%
dplyr::filter(Time == 60) %>%
# Clean up cell value names by removing any strings including and after "__"
cutCellNames(col = "Beta.Metric", sep = "__") %>%
# Filter out raw data, leave only normalized data
dplyr::filter(Data.Type != "Raw") %>%
# Run Tukey test on GLM models
run_tukey_glm_beta(., "Displacement", "Beta.Metric", c("Treatment"),
family_str = "gaussian",
group_by_var = NULL)
| GLM Results - 60 DPE | |||||
| glm(Displacement ~ Treatment); Parasite exposure only at 60 DPE | |||||
| term | estimate | std.error | statistic | p.value | p.adj.sig |
|---|---|---|---|---|---|
| Bray-Curtis_norm | |||||
| (Intercept) | −0.060 | 0.036 | −1.654 | 0.103 | ns |
| A- T- P+ | 0.284 | 0.051 | 5.548 | <0.001 | **** |
| Canberra_norm | |||||
| (Intercept) | −0.124 | 0.028 | −4.457 | <0.001 | **** |
| A- T- P+ | 0.134 | 0.039 | 3.417 | 0.001 | ** |
| ANOVA of GLM - 60 DPE | ||||
| ANOVA(GLM(Displacement ~ Treatment), type = 2); Parasite exposure only at 60 DPE | ||||
| term | statistic | df | p.value | sig |
|---|---|---|---|---|
| Bray-Curtis_norm | ||||
| Treatment | 30.779 | 1.000 | <0.001 | **** |
| Canberra_norm | ||||
| Treatment | 11.677 | 1.000 | <0.001 | *** |
| Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE | ||||||||
| Tukey(Displacement ~ Treatment); Parasite exposure only at 60 DPE | ||||||||
| term | .y. | group1 | group2 | estimate | std.error | statistic | adj.p.value | Variable |
|---|---|---|---|---|---|---|---|---|
| Bray-Curtis_norm | ||||||||
| Treatment | Beta.Score | A- T- P+ | A- T- P- | 0.284 | 0.051 | 5.548 | <0.001 | Treatment |
| Canberra_norm | ||||||||
| Treatment | Beta.Score | A- T- P+ | A- T- P- | 0.134 | 0.039 | 3.417 | <0.001 | Treatment |
| Significant Results:Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE | ||||||||
| Tukey(Displacement ~ Treatment); All treatments at 60 DPE | ||||||||
| term | .y. | group1 | group2 | estimate | std.error | statistic | adj.p.value | Variable |
|---|---|---|---|---|---|---|---|---|
| Bray-Curtis_norm | ||||||||
| Treatment | Beta.Score | A- T- P+ | A- T- P- | 0.284 | 0.051 | 5.548 | <0.001 | Treatment |
| Canberra_norm | ||||||||
| Treatment | Beta.Score | A- T- P+ | A- T- P- | 0.134 | 0.039 | 3.417 | <0.001 | Treatment |
# Set the subsection for this analysis
tmp.resSubSection <- "ParasiteExp_60DPE"
#### Dispersion --------------------------------------------------------------
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]] <-
run_BetaDispersion(dist.matrix = beta.dist.mat[[tmp.resSubSection]],
var = c("treatment_group"))
##### ANOVA ------------------------------------------------------------------
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA"]] <-
get_HoD_anova(betaDisper = beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]],
var = c("treatment_group"))
### Table
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA.Table"]] <-
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA"]] %>%
# Create the table
dplyr::group_by(Beta.Metric) %>%
set_GT(var = "p.value", group.by = "Beta.Metric") %>%
# Title/caption
gt::tab_header(
title = "ANOVA: Homogeneity of Dispersion",
subtitle = "ANOVA(Beta Disperson ~ Treatment); Parasite exposure only at 60 DPE"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "lightgreen"),
locations = gt::cells_body(
columns = c("p.value", "sig"),
rows = p.value < 0.05
)
)
##### Tukey ------------------------------------------------------------------
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey"]] <-
get_HoD_tukey(betaDisper = beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]],
var = c("treatment_group"))
### Table
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey.Table"]] <-
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey"]] %>%
dplyr::mutate(
group1 = case_when(
group1 == "Control" ~ "A- T- P-",
group1 == "Parasite" ~ "A- T- P+",
group1 == "Antibiotics" ~ "A+ T- P-",
group1 == "Antibiotics_Parasite" ~ "A+ T- P+",
group1 == "Temperature" ~ "A- T+ P-",
group1 == "Temperature_Parasite" ~ "A- T+ P+",
group1 == "Antibiotics_Temperature" ~ "A+ T+ P-",
group1 == "Antibiotics_Temperature_Parasite" ~ "A+ T+ P+",
),
group2 = case_when(
group2 == "Control" ~ "A- T- P-",
group2 == "Parasite" ~ "A- T- P+",
group2 == "Antibiotics" ~ "A+ T- P-",
group2 == "Antibiotics_Parasite" ~ "A+ T- P+",
group2 == "Temperature" ~ "A- T+ P-",
group2 == "Temperature_Parasite" ~ "A- T+ P+",
group2 == "Antibiotics_Temperature" ~ "A+ T+ P-",
group2 == "Antibiotics_Temperature_Parasite" ~ "A+ T+ P+",
)) %>%
# Create the table
dplyr::group_by(Beta.Metric) %>%
set_GT(var = "adj.p.value", group.by = "Beta.Metric") %>%
# Title/caption
gt::tab_header(
title = "Tukey: Homogeneity of Dispersion",
subtitle = "Tukey(Beta Disperson ~ Treatment); Parasite exposure only 60 DPE"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "lightgreen"),
locations = gt::cells_body(
columns = c("adj.p.value", "sig"),
rows = adj.p.value < 0.05
)
)
| ANOVA: Homogeneity of Dispersion | ||||||
| ANOVA(Beta Disperson ~ Treatment); Parasite exposure only at 60 DPE | ||||||
| term | df | sumsq | meansq | statistic | p.value | sig |
|---|---|---|---|---|---|---|
| bray | ||||||
| treatment_group | 1.000 | 0.011 | 0.011 | 0.388 | ≥0.25 | ns |
| Residual | 71.000 | 2.088 | 0.029 | NA | NA | NA |
| canberra | ||||||
| treatment_group | 1.000 | 0.001 | 0.001 | 0.608 | ≥0.25 | ns |
| Residual | 71.000 | 0.150 | 0.002 | NA | NA | NA |
| Tukey: Homogeneity of Dispersion | ||||||||
| Tukey(Beta Disperson ~ Treatment); Parasite exposure only 60 DPE | ||||||||
| .y. | term | group1 | group2 | estimate | conf.low | conf.high | adj.p.value | sig |
|---|---|---|---|---|---|---|---|---|
| bray | ||||||||
| Distance | treatment_group | A- T- P+ | A- T- P- | −0.025 | −0.105 | 0.055 | ≥0.25 | ns |
| canberra | ||||||||
| Distance | treatment_group | A- T- P+ | A- T- P- | 0.008 | −0.013 | 0.030 | ≥0.25 | ns |
# Set the subsection for this analysis
tmp.resSubSection <- "ParasiteExp_60DPE"
#### Capscale ----------------------------------------------------------------
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.mod"]] <-
run_capscale(ps.cleaned_alpha %>% # is the phyloseq object we are looping over
microViz::ps_filter(Treatment %in% c(
"A- T- P-",
"A- T- P+"#,
# "A+ T- P-",
# "A+ T- P+",
# "A- T+ P-",
# "A- T+ P+",
# "A+ T+ P-",
# "A+ T+ P+"
)) %>%
microViz::ps_filter(Time == 60),
dist.matrix = beta.dist.mat[[tmp.resSubSection]],
formula_str = "dist ~ Treatment")
##### ADONIS ------------------------------------------------------------------
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS"]] <-
run_cap_adonis(ps.cleaned_alpha %>% # is the phyloseq object we are looping over
microViz::ps_filter(Treatment %in% c(
"A- T- P-",
"A- T- P+"#,
# "A+ T- P-",
# "A+ T- P+",
# "A- T+ P-",
# "A- T+ P+",
# "A+ T+ P-",
# "A+ T+ P+"
)) %>%
microViz::ps_filter(Time == 60),
dist.matrix = beta.dist.mat[[tmp.resSubSection]],
formula_str = "dist ~ Treatment",
by.method = "terms")
## Warning in tidy.anova(.): The following column names in ANOVA output were not
## recognized or transformed: SumOfSqs, R2
## Warning in tidy.anova(.): The following column names in ANOVA output were not
## recognized or transformed: SumOfSqs, R2
### Table
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS.Table"]] <-
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS"]] %>%
set_GT(var = "p.value", group.by = "Beta.Metric") %>%
# Title/caption
gt::tab_header(
title = "ADONIS2",
subtitle = "adonis2(Beta Distance ~ Treatment); Parasite exposure only at 60 DPE"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "lightgreen"),
locations = gt::cells_body(
columns = c("p.value", "sig"),
rows = p.value < 0.05
)
)
| ADONIS2 | ||||||
| adonis2(Beta Distance ~ Treatment); Parasite exposure only at 60 DPE | ||||||
| term | df | SumOfSqs | R2 | statistic | p.value | sig |
|---|---|---|---|---|---|---|
| bray | ||||||
| Treatment | 1.000 | 0.587 | 0.047 | 3.468 | 0.007 | ** |
| Residual | 71.000 | 12.011 | 0.953 | NA | NA | NA |
| Total | 72.000 | 12.597 | 1.000 | NA | NA | NA |
| canberra | ||||||
| Treatment | 1.000 | 1.287 | 0.050 | 3.741 | 0.001 | *** |
| Residual | 71.000 | 24.425 | 0.950 | NA | NA | NA |
| Total | 72.000 | 25.712 | 1.000 | NA | NA | NA |
# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressParaExp"
# GLM Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]] <-
test.alpha.all_metrics$volatility %>%
# Filter
dplyr::filter(Treatment %in% c(
# "A- T- P-",
"A- T- P+",
# "A+ T- P-",
"A+ T- P+",
# "A- T+ P-",
"A- T+ P+",
# "A+ T+ P-",
"A+ T+ P+"
)) %>%
# Clean up cell value names by removing any strings including and after "__"
cutCellNames(col = "Alpha.Metric", sep = "__") %>%
# Filter out raw data, leave only normalized data
dplyr::filter(Data.Type != "Raw") %>%
# Run GLM models
run_glm_models(formula_str = "Displacement ~ Treatment",
alpha_metric_col = "Alpha.Metric",
alpha_score_col = "Displacement",
family_str = "gaussian")
# ANOVA Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["ANOVA"]] <-
run_glm_anova(alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]])
# Tukey Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["Tukey"]] <-
test.alpha.all_metrics$volatility %>%
# Filter
dplyr::filter(Treatment %in% c(
# "A- T- P-",
"A- T- P+",
# "A+ T- P-",
"A+ T- P+",
# "A- T+ P-",
"A- T+ P+",
# "A+ T+ P-",
"A+ T+ P+"
)) %>%
# Clean up cell value names by removing any strings including and after "__"
cutCellNames(col = "Alpha.Metric", sep = "__") %>%
# Filter out raw data, leave only normalized data
dplyr::filter(Data.Type != "Raw") %>%
# Run Tukey test on GLM models
run_tukey_glm(., "Displacement", "Alpha.Metric", c("Treatment"),
family_str = "gaussian",
group_by_var = NULL)
| GLM Results | |||||
| glm(Displacement ~ Treatment); Prior stressors and parasite exposure | |||||
| term | estimate | std.error | statistic | p.value | p.adj.sig |
|---|---|---|---|---|---|
| Shannon | |||||
| (Intercept) | 0.033 | 0.020 | 1.674 | 0.095 | ns |
| A+ T- P+ | −0.027 | 0.028 | −0.956 | ≥0.25 | ns |
| A- T+ P+ | −0.047 | 0.028 | −1.675 | 0.095 | ns |
| A+ T+ P+ | −0.017 | 0.029 | −0.575 | ≥0.25 | ns |
| Simpson | |||||
| (Intercept) | −0.006 | 0.028 | −0.219 | ≥0.25 | ns |
| A+ T- P+ | −0.037 | 0.040 | −0.916 | ≥0.25 | ns |
| A- T+ P+ | −0.105 | 0.041 | −2.592 | 0.010 | * |
| A+ T+ P+ | −0.042 | 0.042 | −1.020 | ≥0.25 | ns |
| Richness | |||||
| (Intercept) | 0.230 | 0.025 | 9.206 | <0.001 | **** |
| A+ T- P+ | −0.054 | 0.036 | −1.519 | 0.130 | ns |
| A- T+ P+ | 0.007 | 0.036 | 0.186 | ≥0.25 | ns |
| A+ T+ P+ | −0.063 | 0.037 | −1.711 | 0.088 | ns |
| ANOVA of GLM | ||||
| ANOVA(GLM(Displacement ~ Treatment), type = 2); Prior stressors and parasite exposure | ||||
| term | statistic | df | p.value | sig |
|---|---|---|---|---|
| Shannon | ||||
| Treatment | 2.931 | 3.000 | ≥0.25 | ns |
| Simpson | ||||
| Treatment | 6.896 | 3.000 | 0.075 | ns |
| Richness | ||||
| Treatment | 5.804 | 3.000 | 0.122 | ns |
| Pairwise Tukey's HSD, p.adj: Dunnett | ||||||||
| Tukey(Displacement ~ Treatment); Prior stressors and parasite exposure | ||||||||
| term | .y. | group1 | group2 | estimate | std.error | statistic | adj.p.value | Variable |
|---|---|---|---|---|---|---|---|---|
| Shannon | ||||||||
| Treatment | Alpha.Score | A+ T- P+ | A- T- P+ | −0.027 | 0.028 | −0.956 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A- T- P+ | −0.047 | 0.028 | −1.675 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T- P+ | −0.017 | 0.029 | −0.575 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A+ T- P+ | −0.020 | 0.028 | −0.721 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A+ T- P+ | 0.010 | 0.029 | 0.349 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T+ P+ | 0.031 | 0.029 | 1.046 | ≥0.25 | Treatment |
| Simpson | ||||||||
| Treatment | Alpha.Score | A+ T- P+ | A- T- P+ | −0.037 | 0.040 | −0.916 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A- T- P+ | −0.105 | 0.041 | −2.592 | 0.047 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T- P+ | −0.042 | 0.042 | −1.020 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A+ T- P+ | −0.068 | 0.041 | −1.671 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A+ T- P+ | −0.005 | 0.042 | −0.131 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T+ P+ | 0.063 | 0.042 | 1.490 | ≥0.25 | Treatment |
| Richness | ||||||||
| Treatment | Alpha.Score | A+ T- P+ | A- T- P+ | −0.054 | 0.036 | −1.519 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A- T- P+ | 0.007 | 0.036 | 0.186 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T- P+ | −0.063 | 0.037 | −1.711 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A+ T- P+ | 0.061 | 0.036 | 1.685 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A+ T- P+ | −0.009 | 0.037 | −0.237 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T+ P+ | −0.070 | 0.037 | −1.870 | 0.241 | Treatment |
# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressParaExp_60DPE"
# GLM Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]] <-
test.alpha.all_metrics$displacement %>%
dplyr::filter(Treatment %in% c(
# "A- T- P-",
"A- T- P+",
# "A+ T- P-",
"A+ T- P+",
# "A- T+ P-",
"A- T+ P+",
# "A+ T+ P-",
"A+ T+ P+"
)) %>%
dplyr::filter(Time == 60) %>%
# Clean up cell value names by removing any strings including and after "__"
cutCellNames(col = "Alpha.Metric", sep = "__") %>%
# Filter out raw data, leave only normalized data
dplyr::filter(Data.Type != "Raw") %>%
# Run GLM models
run_glm_models(formula_str = "Displacement ~ Treatment",
alpha_metric_col = "Alpha.Metric",
alpha_score_col = "Displacement",
family_str = "gaussian")
# ANOVA Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["ANOVA"]] <-
run_glm_anova(alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]])
# Tukey Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["Tukey"]] <-
test.alpha.all_metrics$displacement %>%
dplyr::filter(Treatment %in% c(
# "A- T- P-",
"A- T- P+",
# "A+ T- P-",
"A+ T- P+",
# "A- T+ P-",
"A- T+ P+",
# "A+ T+ P-",
"A+ T+ P+"
)) %>%
dplyr::filter(Time == 60) %>%
# Clean up cell value names by removing any strings including and after "__"
cutCellNames(col = "Alpha.Metric", sep = "__") %>%
# Filter out raw data, leave only normalized data
dplyr::filter(Data.Type != "Raw") %>%
# Run Tukey test on GLM models
run_tukey_glm(., "Displacement", "Alpha.Metric", c("Treatment"),
family_str = "gaussian",
group_by_var = NULL)
| GLM Results - 60 DPE | |||||
| glm(Displacement ~ Treatment); All treatments at 60 DPE | |||||
| term | estimate | std.error | statistic | p.value | p.adj.sig |
|---|---|---|---|---|---|
| Shannon | |||||
| (Intercept) | 0.008 | 0.037 | 0.224 | ≥0.25 | ns |
| A+ T- P+ | −0.001 | 0.053 | −0.022 | ≥0.25 | ns |
| A- T+ P+ | −0.077 | 0.054 | −1.424 | 0.157 | ns |
| A+ T+ P+ | 0.004 | 0.058 | 0.065 | ≥0.25 | ns |
| Simpson | |||||
| (Intercept) | −0.062 | 0.048 | −1.286 | 0.201 | ns |
| A+ T- P+ | 0.000 | 0.069 | −0.005 | ≥0.25 | ns |
| A- T+ P+ | −0.180 | 0.071 | −2.547 | 0.012 | * |
| A+ T+ P+ | −0.037 | 0.076 | −0.484 | ≥0.25 | ns |
| Richness | |||||
| (Intercept) | 0.363 | 0.032 | 11.248 | <0.001 | **** |
| A+ T- P+ | −0.132 | 0.046 | −2.863 | 0.005 | ** |
| A- T+ P+ | −0.003 | 0.047 | −0.069 | ≥0.25 | ns |
| A+ T+ P+ | −0.053 | 0.051 | −1.039 | ≥0.25 | ns |
| ANOVA of GLM - 60 DPE | ||||
| ANOVA(GLM(Displacement ~ Treatment), type = 2); Prior stressors and parasite exposure at 60 DPE | ||||
| term | statistic | df | p.value | sig |
|---|---|---|---|---|
| Shannon | ||||
| Treatment | 2.892 | 3.000 | ≥0.25 | ns |
| Simpson | ||||
| Treatment | 8.445 | 3.000 | 0.038 | * |
| Richness | ||||
| Treatment | 10.375 | 3.000 | 0.016 | * |
| Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE | ||||||||
| Tukey(Displacement ~ Treatment); Prior stressors and parasite exposure at 60 DPE | ||||||||
| term | .y. | group1 | group2 | estimate | std.error | statistic | adj.p.value | Variable |
|---|---|---|---|---|---|---|---|---|
| Shannon | ||||||||
| Treatment | Alpha.Score | A+ T- P+ | A- T- P+ | −0.001 | 0.053 | −0.022 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A- T- P+ | −0.077 | 0.054 | −1.424 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T- P+ | 0.004 | 0.058 | 0.065 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A+ T- P+ | −0.076 | 0.055 | −1.384 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A+ T- P+ | 0.005 | 0.059 | 0.084 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T+ P+ | 0.081 | 0.060 | 1.351 | ≥0.25 | Treatment |
| Simpson | ||||||||
| Treatment | Alpha.Score | A+ T- P+ | A- T- P+ | 0.000 | 0.069 | −0.005 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A- T- P+ | −0.180 | 0.071 | −2.547 | 0.053 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T- P+ | −0.037 | 0.076 | −0.484 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A+ T- P+ | −0.179 | 0.071 | −2.509 | 0.058 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A+ T- P+ | −0.036 | 0.077 | −0.475 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T+ P+ | 0.143 | 0.078 | 1.834 | ≥0.25 | Treatment |
| Richness | ||||||||
| Treatment | Alpha.Score | A+ T- P+ | A- T- P+ | −0.132 | 0.046 | −2.863 | 0.022 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A- T- P+ | −0.003 | 0.047 | −0.069 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T- P+ | −0.053 | 0.051 | −1.039 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A+ T- P+ | 0.129 | 0.048 | 2.692 | 0.036 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A+ T- P+ | 0.080 | 0.051 | 1.551 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P+ | A- T+ P+ | −0.050 | 0.052 | −0.946 | ≥0.25 | Treatment |
| Significant Results:Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE | ||||||||
| Tukey(Displacement ~ Treatment); All treatments at 60 DPE | ||||||||
| term | .y. | group1 | group2 | estimate | std.error | statistic | adj.p.value | Variable |
|---|---|---|---|---|---|---|---|---|
| Richness | ||||||||
| Treatment | Alpha.Score | A+ T- P+ | A- T- P+ | −0.132 | 0.046 | −2.863 | 0.022 | Treatment |
| Treatment | Alpha.Score | A- T+ P+ | A+ T- P+ | 0.129 | 0.048 | 2.692 | 0.036 | Treatment |
# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressParaExp"
# GLM Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]] <-
test.beta.all_metrics$volatility %>%
# Filter
dplyr::filter(Treatment %in% c(
# "A- T- P-",
"A- T- P+",
# "A+ T- P-",
"A+ T- P+",
# "A- T+ P-",
"A- T+ P+",
# "A+ T+ P-",
"A+ T+ P+"
)) %>%
# Clean up cell value names by removing any strings including and after "__"
cutCellNames(col = "Beta.Metric", sep = "__") %>%
# Filter out raw data, leave only normalized data
dplyr::filter(Data.Type != "Raw") %>%
# Run GLM models
run_glm_models_beta(formula_str = "Displacement ~ Treatment",
beta_metric_col = "Beta.Metric",
beta_score_col = "Displacement",
family_str = "gaussian")
# ANOVA Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["ANOVA"]] <-
run_glm_anova_beta(beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]])
# Tukey Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["Tukey"]] <-
test.beta.all_metrics$volatility %>%
# Filter
dplyr::filter(Treatment %in% c(
# "A- T- P-",
"A- T- P+",
# "A+ T- P-",
"A+ T- P+",
# "A- T+ P-",
"A- T+ P+",
# "A+ T+ P-",
"A+ T+ P+"
)) %>%
# Clean up cell value names by removing any strings including and after "__"
cutCellNames(col = "Beta.Metric", sep = "__") %>%
# Filter out raw data, leave only normalized data
dplyr::filter(Data.Type != "Raw") %>%
# Run Tukey test on GLM models
run_tukey_glm_beta(., "Displacement", "Beta.Metric", c("Treatment"),
family_str = "gaussian",
group_by_var = NULL)
| GLM Results - Beta | |||||
| glm(Displacement ~ Treatment); Prior stressors and parasite exposure | |||||
| term | estimate | std.error | statistic | p.value | p.adj.sig |
|---|---|---|---|---|---|
| Bray-Curtis_norm | |||||
| (Intercept) | 0.081 | 0.027 | 2.973 | 0.003 | ** |
| A+ T- P+ | −0.141 | 0.039 | −3.633 | <0.001 | *** |
| A- T+ P+ | −0.306 | 0.039 | −7.843 | <0.001 | **** |
| A+ T+ P+ | −0.146 | 0.040 | −3.662 | <0.001 | *** |
| Canberra_norm | |||||
| (Intercept) | −0.106 | 0.020 | −5.372 | <0.001 | **** |
| A+ T- P+ | 0.013 | 0.028 | 0.478 | ≥0.25 | ns |
| A- T+ P+ | −0.016 | 0.028 | −0.573 | ≥0.25 | ns |
| A+ T+ P+ | 0.026 | 0.029 | 0.900 | ≥0.25 | ns |
| ANOVA of GLM - Beta | ||||
| ANOVA(GLM(Displacement ~ Treatment), type = 2); Prior stressors and parasite exposure | ||||
| term | statistic | df | p.value | sig |
|---|---|---|---|---|
| Bray-Curtis_norm | ||||
| Treatment | 61.580 | 3.000 | <0.001 | **** |
| Canberra_norm | ||||
| Treatment | 2.326 | 3.000 | ≥0.25 | ns |
| Pairwise Tukey's HSD, p.adj: Dunnett - Beta | ||||||||
| Tukey(Displacement ~ Treatment); Prior stressors and parasite exposure | ||||||||
| term | .y. | group1 | group2 | estimate | std.error | statistic | adj.p.value | Variable |
|---|---|---|---|---|---|---|---|---|
| Bray-Curtis_norm | ||||||||
| Treatment | Beta.Score | A+ T- P+ | A- T- P+ | −0.141 | 0.039 | −3.633 | 0.002 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A- T- P+ | −0.306 | 0.039 | −7.843 | <0.001 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A- T- P+ | −0.146 | 0.040 | −3.662 | 0.002 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A+ T- P+ | −0.165 | 0.039 | −4.208 | <0.001 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A+ T- P+ | −0.006 | 0.040 | −0.141 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A- T+ P+ | 0.159 | 0.040 | 3.940 | <0.001 | Treatment |
| Canberra_norm | ||||||||
| Treatment | Beta.Score | A+ T- P+ | A- T- P+ | 0.013 | 0.028 | 0.478 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A- T- P+ | −0.016 | 0.028 | −0.573 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A- T- P+ | 0.026 | 0.029 | 0.900 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A+ T- P+ | −0.030 | 0.029 | −1.042 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A+ T- P+ | 0.013 | 0.029 | 0.434 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A- T+ P+ | 0.042 | 0.029 | 1.442 | ≥0.25 | Treatment |
# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressParaExp_60DPE"
# GLM Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]] <-
test.beta.all_metrics$displacement %>%
dplyr::filter(Treatment %in% c(
# "A- T- P-",
"A- T- P+",
# "A+ T- P-",
"A+ T- P+",
# "A- T+ P-",
"A- T+ P+",
# "A+ T+ P-",
"A+ T+ P+"
)) %>%
dplyr::filter(Time == 60) %>%
# Clean up cell value names by removing any strings including and after "__"
cutCellNames(col = "Beta.Metric", sep = "__") %>%
# Filter out raw data, leave only normalized data
dplyr::filter(Data.Type != "Raw") %>%
# Run GLM models
run_glm_models_beta(formula_str = "Displacement ~ Treatment",
beta_metric_col = "Beta.Metric",
beta_score_col = "Displacement",
family_str = "gaussian")
# ANOVA Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["ANOVA"]] <-
run_glm_anova_beta(beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]])
# Tukey Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["Tukey"]] <-
test.beta.all_metrics$displacement %>%
dplyr::filter(Treatment %in% c(
# "A- T- P-",
"A- T- P+",
# "A+ T- P-",
"A+ T- P+",
# "A- T+ P-",
"A- T+ P+",
# "A+ T+ P-",
"A+ T+ P+"
)) %>%
dplyr::filter(Time == 60) %>%
# Clean up cell value names by removing any strings including and after "__"
cutCellNames(col = "Beta.Metric", sep = "__") %>%
# Filter out raw data, leave only normalized data
dplyr::filter(Data.Type != "Raw") %>%
# Run Tukey test on GLM models
run_tukey_glm_beta(., "Displacement", "Beta.Metric", c("Treatment"),
family_str = "gaussian",
group_by_var = NULL)
| GLM Results - 60 DPE | |||||
| glm(Displacement ~ Treatment); Prior stressors and parasite exposure only at 60 DPE | |||||
| term | estimate | std.error | statistic | p.value | p.adj.sig |
|---|---|---|---|---|---|
| Bray-Curtis_norm | |||||
| (Intercept) | 0.224 | 0.038 | 5.823 | <0.001 | **** |
| A+ T- P+ | −0.305 | 0.055 | −5.534 | <0.001 | **** |
| A- T+ P+ | −0.490 | 0.056 | −8.688 | <0.001 | **** |
| A+ T+ P+ | −0.220 | 0.060 | −3.646 | <0.001 | *** |
| Canberra_norm | |||||
| (Intercept) | 0.010 | 0.028 | 0.340 | ≥0.25 | ns |
| A+ T- P+ | −0.097 | 0.040 | −2.398 | 0.018 | * |
| A- T+ P+ | −0.069 | 0.041 | −1.670 | 0.097 | ns |
| A+ T+ P+ | −0.063 | 0.044 | −1.425 | 0.157 | ns |
| ANOVA of GLM - 60 DPE | ||||
| ANOVA(GLM(Displacement ~ Treatment), type = 2); Prior stressors and parasite exposure only at 60 DPE | ||||
| term | statistic | df | p.value | sig |
|---|---|---|---|---|
| Bray-Curtis_norm | ||||
| Treatment | 78.442 | 3.000 | <0.001 | **** |
| Canberra_norm | ||||
| Treatment | 6.174 | 3.000 | 0.103 | ns |
| Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE | ||||||||
| Tukey(Displacement ~ Treatment); Prior stressors and parasite exposure only at 60 DPE | ||||||||
| term | .y. | group1 | group2 | estimate | std.error | statistic | adj.p.value | Variable |
|---|---|---|---|---|---|---|---|---|
| Bray-Curtis_norm | ||||||||
| Treatment | Beta.Score | A+ T- P+ | A- T- P+ | −0.305 | 0.055 | −5.534 | <0.001 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A- T- P+ | −0.490 | 0.056 | −8.688 | 0.000 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A- T- P+ | −0.220 | 0.060 | −3.646 | 0.002 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A+ T- P+ | −0.185 | 0.057 | −3.240 | 0.007 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A+ T- P+ | 0.084 | 0.061 | 1.379 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A- T+ P+ | 0.269 | 0.062 | 4.321 | <0.001 | Treatment |
| Canberra_norm | ||||||||
| Treatment | Beta.Score | A+ T- P+ | A- T- P+ | −0.097 | 0.040 | −2.398 | 0.077 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A- T- P+ | −0.069 | 0.041 | −1.670 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A- T- P+ | −0.063 | 0.044 | −1.425 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A+ T- P+ | 0.028 | 0.042 | 0.663 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A+ T- P+ | 0.034 | 0.045 | 0.750 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A- T+ P+ | 0.006 | 0.046 | 0.128 | ≥0.25 | Treatment |
| Significant Results:Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE | ||||||||
| Tukey(Displacement ~ Treatment); All treatments at 60 DPE | ||||||||
| term | .y. | group1 | group2 | estimate | std.error | statistic | adj.p.value | Variable |
|---|---|---|---|---|---|---|---|---|
| Bray-Curtis_norm | ||||||||
| Treatment | Beta.Score | A+ T- P+ | A- T- P+ | −0.305 | 0.055 | −5.534 | <0.001 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A- T- P+ | −0.490 | 0.056 | −8.688 | 0.000 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A- T- P+ | −0.220 | 0.060 | −3.646 | 0.002 | Treatment |
| Treatment | Beta.Score | A- T+ P+ | A+ T- P+ | −0.185 | 0.057 | −3.240 | 0.007 | Treatment |
| Treatment | Beta.Score | A+ T+ P+ | A- T+ P+ | 0.269 | 0.062 | 4.321 | <0.001 | Treatment |
# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressParaExp_60DPE"
#### Dispersion --------------------------------------------------------------
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]] <-
run_BetaDispersion(dist.matrix = beta.dist.mat[[tmp.resSubSection]],
var = c("treatment_group"))
##### ANOVA ------------------------------------------------------------------
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA"]] <-
get_HoD_anova(betaDisper = beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]],
var = c("treatment_group"))
### Table
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA.Table"]] <-
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA"]] %>%
# Create the table
dplyr::group_by(Beta.Metric) %>%
set_GT(var = "p.value", group.by = "Beta.Metric") %>%
# Title/caption
gt::tab_header(
title = "ANOVA: Homogeneity of Dispersion",
subtitle = "ANOVA(Beta Disperson ~ Treatment); Prior stressors and parasite exposure at 60 DPE"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "lightgreen"),
locations = gt::cells_body(
columns = c("p.value", "sig"),
rows = p.value < 0.05
)
)
##### Tukey ------------------------------------------------------------------
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey"]] <-
get_HoD_tukey(betaDisper = beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]],
var = c("treatment_group"))
### Table
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey.Table"]] <-
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey"]] %>%
dplyr::mutate(
group1 = case_when(
group1 == "Control" ~ "A- T- P-",
group1 == "Parasite" ~ "A- T- P+",
group1 == "Antibiotics" ~ "A+ T- P-",
group1 == "Antibiotics_Parasite" ~ "A+ T- P+",
group1 == "Temperature" ~ "A- T+ P-",
group1 == "Temperature_Parasite" ~ "A- T+ P+",
group1 == "Antibiotics_Temperature" ~ "A+ T+ P-",
group1 == "Antibiotics_Temperature_Parasite" ~ "A+ T+ P+",
),
group2 = case_when(
group2 == "Control" ~ "A- T- P-",
group2 == "Parasite" ~ "A- T- P+",
group2 == "Antibiotics" ~ "A+ T- P-",
group2 == "Antibiotics_Parasite" ~ "A+ T- P+",
group2 == "Temperature" ~ "A- T+ P-",
group2 == "Temperature_Parasite" ~ "A- T+ P+",
group2 == "Antibiotics_Temperature" ~ "A+ T+ P-",
group2 == "Antibiotics_Temperature_Parasite" ~ "A+ T+ P+",
)) %>%
# Create the table
dplyr::group_by(Beta.Metric) %>%
set_GT(var = "adj.p.value", group.by = "Beta.Metric") %>%
# Title/caption
gt::tab_header(
title = "Tukey: Homogeneity of Dispersion",
subtitle = "Tukey(Beta Disperson ~ Treatment); Prior stressors and parasite exposure at 60 DPE"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "lightgreen"),
locations = gt::cells_body(
columns = c("adj.p.value", "sig"),
rows = adj.p.value < 0.05
)
)
| ANOVA: Homogeneity of Dispersion | ||||||
| ANOVA(Beta Disperson ~ Treatment); Prior stressors and parasite exposure at 60 DPE | ||||||
| term | df | sumsq | meansq | statistic | p.value | sig |
|---|---|---|---|---|---|---|
| bray | ||||||
| treatment_group | 3.000 | 0.291 | 0.097 | 3.785 | 0.012 | * |
| Residual | 125.000 | 3.198 | 0.026 | NA | NA | NA |
| canberra | ||||||
| treatment_group | 3.000 | 0.049 | 0.016 | 4.464 | 0.005 | ** |
| Residual | 125.000 | 0.462 | 0.004 | NA | NA | NA |
| Tukey: Homogeneity of Dispersion | ||||||||
| Tukey(Beta Disperson ~ Treatment); Prior stressors and parasite exposure at 60 DPE | ||||||||
| .y. | term | group1 | group2 | estimate | conf.low | conf.high | adj.p.value | sig |
|---|---|---|---|---|---|---|---|---|
| bray | ||||||||
| Distance | treatment_group | A+ T- P+ | A- T- P+ | −0.040 | −0.138 | 0.058 | ≥0.25 | ns |
| Distance | treatment_group | A- T+ P+ | A- T- P+ | −0.081 | −0.182 | 0.019 | 0.157 | ns |
| Distance | treatment_group | A+ T+ P+ | A- T- P+ | −0.132 | −0.240 | −0.024 | 0.009 | ** |
| Distance | treatment_group | A- T+ P+ | A+ T- P+ | −0.041 | −0.143 | 0.061 | ≥0.25 | ns |
| Distance | treatment_group | A+ T+ P+ | A+ T- P+ | −0.092 | −0.201 | 0.017 | 0.130 | ns |
| Distance | treatment_group | A+ T+ P+ | A- T+ P+ | −0.051 | −0.162 | 0.060 | ≥0.25 | ns |
| canberra | ||||||||
| Distance | treatment_group | A+ T- P+ | A- T- P+ | −0.001 | −0.038 | 0.036 | ≥0.25 | ns |
| Distance | treatment_group | A- T+ P+ | A- T- P+ | −0.015 | −0.053 | 0.023 | ≥0.25 | ns |
| Distance | treatment_group | A+ T+ P+ | A- T- P+ | −0.052 | −0.093 | −0.011 | 0.007 | ** |
| Distance | treatment_group | A- T+ P+ | A+ T- P+ | −0.014 | −0.053 | 0.025 | ≥0.25 | ns |
| Distance | treatment_group | A+ T+ P+ | A+ T- P+ | −0.051 | −0.093 | −0.010 | 0.009 | ** |
| Distance | treatment_group | A+ T+ P+ | A- T+ P+ | −0.037 | −0.079 | 0.005 | 0.106 | ns |
# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressParaExp_60DPE"
#### Capscale ----------------------------------------------------------------
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.mod"]] <-
run_capscale(ps.cleaned_alpha %>% # is the phyloseq object we are looping over
microViz::ps_filter(Treatment %in% c(
# "A- T- P-",
"A- T- P+",
# "A+ T- P-",
"A+ T- P+",
# "A- T+ P-",
"A- T+ P+",
# "A+ T+ P-",
"A+ T+ P+"
)) %>%
microViz::ps_filter(Time == 60),
dist.matrix = beta.dist.mat[[tmp.resSubSection]],
formula_str = "dist ~ Treatment")
##### ADONIS ------------------------------------------------------------------
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS"]] <-
run_cap_adonis(ps.cleaned_alpha %>% # is the phyloseq object we are looping over
microViz::ps_filter(Treatment %in% c(
# "A- T- P-",
"A- T- P+",
# "A+ T- P-",
"A+ T- P+",
# "A- T+ P-",
"A- T+ P+",
# "A+ T+ P-",
"A+ T+ P+"
)) %>%
microViz::ps_filter(Time == 60),
dist.matrix = beta.dist.mat[[tmp.resSubSection]],
formula_str = "dist ~ Treatment",
by.method = "terms")
## Warning in tidy.anova(.): The following column names in ANOVA output were not
## recognized or transformed: SumOfSqs, R2
## Warning in tidy.anova(.): The following column names in ANOVA output were not
## recognized or transformed: SumOfSqs, R2
### Table
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS.Table"]] <-
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS"]] %>%
set_GT(var = "p.value", group.by = "Beta.Metric") %>%
# Title/caption
gt::tab_header(
title = "ADONIS2",
subtitle = "adonis2(Beta Distance ~ Treatment); Prior stressors and parasite exposure at 60 DPE"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "lightgreen"),
locations = gt::cells_body(
columns = c("p.value", "sig"),
rows = p.value < 0.05
)
)
| ADONIS2 | ||||||
| adonis2(Beta Distance ~ Treatment); Prior stressors and parasite exposure at 60 DPE | ||||||
| term | df | SumOfSqs | R2 | statistic | p.value | sig |
|---|---|---|---|---|---|---|
| bray | ||||||
| Treatment | 3.000 | 0.708 | 0.045 | 1.959 | 0.037 | * |
| Residual | 125.000 | 15.053 | 0.955 | NA | NA | NA |
| Total | 128.000 | 15.761 | 1.000 | NA | NA | NA |
| canberra | ||||||
| Treatment | 3.000 | 2.852 | 0.064 | 2.833 | 0.001 | *** |
| Residual | 125.000 | 41.942 | 0.936 | NA | NA | NA |
| Total | 128.000 | 44.794 | 1.000 | NA | NA | NA |
# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressNoParaExp"
# GLM Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]] <-
test.alpha.all_metrics$volatility %>%
# Filter
dplyr::filter(Treatment %in% c(
"A- T- P-",
# "A- T- P+",
"A+ T- P-",
# "A+ T- P+",
"A- T+ P-",
# "A- T+ P+",
"A+ T+ P-"#,
# "A+ T+ P+"
)) %>%
# Clean up cell value names by removing any strings including and after "__"
cutCellNames(col = "Alpha.Metric", sep = "__") %>%
# Filter out raw data, leave only normalized data
dplyr::filter(Data.Type != "Raw") %>%
# Run GLM models
run_glm_models(formula_str = "Displacement ~ Treatment",
alpha_metric_col = "Alpha.Metric",
alpha_score_col = "Displacement",
family_str = "gaussian")
# ANOVA Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["ANOVA"]] <-
run_glm_anova(alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]])
# Tukey Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["Tukey"]] <-
test.alpha.all_metrics$volatility %>%
# Filter
dplyr::filter(Treatment %in% c(
"A- T- P-",
# "A- T- P+",
"A+ T- P-",
# "A+ T- P+",
"A- T+ P-",
# "A- T+ P+",
"A+ T+ P-"#,
# "A+ T+ P+"
)) %>%
# Clean up cell value names by removing any strings including and after "__"
cutCellNames(col = "Alpha.Metric", sep = "__") %>%
# Filter out raw data, leave only normalized data
dplyr::filter(Data.Type != "Raw") %>%
# Run Tukey test on GLM models
run_tukey_glm(., "Displacement", "Alpha.Metric", c("Treatment"),
family_str = "gaussian",
group_by_var = NULL)
| GLM Results | |||||
| glm(Displacement ~ Treatment); Prior stressors and no parasite exposure | |||||
| term | estimate | std.error | statistic | p.value | p.adj.sig |
|---|---|---|---|---|---|
| Shannon | |||||
| (Intercept) | 0.078 | 0.019 | 4.047 | <0.001 | **** |
| A+ T- P- | −0.051 | 0.028 | −1.839 | 0.067 | ns |
| A- T+ P- | −0.004 | 0.028 | −0.153 | ≥0.25 | ns |
| A+ T+ P- | −0.113 | 0.029 | −3.849 | <0.001 | *** |
| Simpson | |||||
| (Intercept) | 0.042 | 0.027 | 1.562 | 0.119 | ns |
| A+ T- P- | −0.034 | 0.039 | −0.861 | ≥0.25 | ns |
| A- T+ P- | 0.035 | 0.040 | 0.894 | ≥0.25 | ns |
| A+ T+ P- | −0.181 | 0.042 | −4.369 | <0.001 | **** |
| Richness | |||||
| (Intercept) | 0.221 | 0.023 | 9.748 | <0.001 | **** |
| A+ T- P- | −0.049 | 0.033 | −1.508 | 0.133 | ns |
| A- T+ P- | −0.016 | 0.033 | −0.493 | ≥0.25 | ns |
| A+ T+ P- | 0.029 | 0.035 | 0.825 | ≥0.25 | ns |
| ANOVA of GLM | ||||
| ANOVA(GLM(Displacement ~ Treatment), type = 2); Prior stressors and no parasite exposure | ||||
| term | statistic | df | p.value | sig |
|---|---|---|---|---|
| Shannon | ||||
| Treatment | 18.423 | 3.000 | <0.001 | *** |
| Simpson | ||||
| Treatment | 29.201 | 3.000 | <0.001 | **** |
| Richness | ||||
| Treatment | 5.193 | 3.000 | 0.158 | ns |
| Pairwise Tukey's HSD, p.adj: Dunnett | ||||||||
| Tukey(Displacement ~ Treatment); Prior stressors and no parasite exposure | ||||||||
| term | .y. | group1 | group2 | estimate | std.error | statistic | adj.p.value | Variable |
|---|---|---|---|---|---|---|---|---|
| Shannon | ||||||||
| Treatment | Alpha.Score | A+ T- P- | A- T- P- | −0.051 | 0.028 | −1.839 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A- T- P- | −0.004 | 0.028 | −0.153 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T- P- | −0.113 | 0.029 | −3.849 | <0.001 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A+ T- P- | 0.047 | 0.029 | 1.629 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A+ T- P- | −0.062 | 0.030 | −2.083 | 0.158 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T+ P- | −0.109 | 0.030 | −3.596 | 0.002 | Treatment |
| Simpson | ||||||||
| Treatment | Alpha.Score | A+ T- P- | A- T- P- | −0.034 | 0.039 | −0.861 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A- T- P- | 0.035 | 0.040 | 0.894 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T- P- | −0.181 | 0.042 | −4.369 | <0.001 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A+ T- P- | 0.069 | 0.040 | 1.711 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A+ T- P- | −0.148 | 0.042 | −3.500 | 0.003 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T+ P- | −0.217 | 0.043 | −5.073 | <0.001 | Treatment |
| Richness | ||||||||
| Treatment | Alpha.Score | A+ T- P- | A- T- P- | −0.049 | 0.033 | −1.508 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A- T- P- | −0.016 | 0.033 | −0.493 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T- P- | 0.029 | 0.035 | 0.825 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A+ T- P- | 0.033 | 0.034 | 0.975 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A+ T- P- | 0.078 | 0.035 | 2.206 | 0.121 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T+ P- | 0.045 | 0.036 | 1.258 | ≥0.25 | Treatment |
# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressNoParaExp_60DPE"
# GLM Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]] <-
test.alpha.all_metrics$displacement %>%
dplyr::filter(Treatment %in% c(
"A- T- P-",
# "A- T- P+",
"A+ T- P-",
# "A+ T- P+",
"A- T+ P-",
# "A- T+ P+",
"A+ T+ P-"#,
# "A+ T+ P+"
)) %>%
dplyr::filter(Time == 60) %>%
# Clean up cell value names by removing any strings including and after "__"
cutCellNames(col = "Alpha.Metric", sep = "__") %>%
# Filter out raw data, leave only normalized data
dplyr::filter(Data.Type != "Raw") %>%
# Run GLM models
run_glm_models(formula_str = "Displacement ~ Treatment",
alpha_metric_col = "Alpha.Metric",
alpha_score_col = "Displacement",
family_str = "gaussian")
# ANOVA Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["ANOVA"]] <-
run_glm_anova(alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]])
# Tukey Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["Tukey"]] <-
test.alpha.all_metrics$displacement %>%
dplyr::filter(Treatment %in% c(
"A- T- P-",
# "A- T- P+",
"A+ T- P-",
# "A+ T- P+",
"A- T+ P-",
# "A- T+ P+",
"A+ T+ P-"#,
# "A+ T+ P+"
)) %>%
dplyr::filter(Time == 60) %>%
# Clean up cell value names by removing any strings including and after "__"
cutCellNames(col = "Alpha.Metric", sep = "__") %>%
# Filter out raw data, leave only normalized data
dplyr::filter(Data.Type != "Raw") %>%
# Run Tukey test on GLM models
run_tukey_glm(., "Displacement", "Alpha.Metric", c("Treatment"),
family_str = "gaussian",
group_by_var = NULL)
| GLM Results - 60 DPE | |||||
| glm(Displacement ~ Treatment); Prior stressors and no parasite exposure at 60 DPE | |||||
| term | estimate | std.error | statistic | p.value | p.adj.sig |
|---|---|---|---|---|---|
| Shannon | |||||
| (Intercept) | 0.091 | 0.036 | 2.496 | 0.014 | * |
| A+ T- P- | −0.034 | 0.054 | −0.627 | ≥0.25 | ns |
| A- T+ P- | −0.041 | 0.056 | −0.731 | ≥0.25 | ns |
| A+ T+ P- | −0.191 | 0.067 | −2.860 | 0.005 | ** |
| Simpson | |||||
| (Intercept) | 0.008 | 0.050 | 0.161 | ≥0.25 | ns |
| A+ T- P- | 0.023 | 0.074 | 0.308 | ≥0.25 | ns |
| A- T+ P- | 0.030 | 0.077 | 0.395 | ≥0.25 | ns |
| A+ T+ P- | −0.332 | 0.092 | −3.627 | <0.001 | *** |
| Richness | |||||
| (Intercept) | 0.338 | 0.031 | 10.867 | <0.001 | **** |
| A+ T- P- | −0.064 | 0.046 | −1.388 | 0.168 | ns |
| A- T+ P- | −0.089 | 0.048 | −1.856 | 0.066 | ns |
| A+ T+ P- | 0.067 | 0.057 | 1.174 | 0.243 | ns |
| ANOVA of GLM - 60 DPE | ||||
| ANOVA(GLM(Displacement ~ Treatment), type = 2); Prior stressors and no parasite exposure at 60 DPE | ||||
| term | statistic | df | p.value | sig |
|---|---|---|---|---|
| Shannon | ||||
| Treatment | 8.403 | 3.000 | 0.038 | * |
| Simpson | ||||
| Treatment | 17.775 | 3.000 | <0.001 | *** |
| Richness | ||||
| Treatment | 8.617 | 3.000 | 0.035 | * |
| Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE | ||||||||
| Tukey(Displacement ~ Treatment); Prior stressors and no parasite exposure at 60 DPE | ||||||||
| term | .y. | group1 | group2 | estimate | std.error | statistic | adj.p.value | Variable |
|---|---|---|---|---|---|---|---|---|
| Shannon | ||||||||
| Treatment | Alpha.Score | A+ T- P- | A- T- P- | −0.034 | 0.054 | −0.627 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A- T- P- | −0.041 | 0.056 | −0.731 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T- P- | −0.191 | 0.067 | −2.860 | 0.022 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A+ T- P- | −0.007 | 0.058 | −0.124 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A+ T- P- | −0.158 | 0.069 | −2.289 | 0.099 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T+ P- | −0.150 | 0.071 | −2.131 | 0.142 | Treatment |
| Simpson | ||||||||
| Treatment | Alpha.Score | A+ T- P- | A- T- P- | 0.023 | 0.074 | 0.308 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A- T- P- | 0.030 | 0.077 | 0.395 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T- P- | −0.332 | 0.092 | −3.627 | 0.002 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A+ T- P- | 0.008 | 0.080 | 0.095 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A+ T- P- | −0.355 | 0.094 | −3.765 | <0.001 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T+ P- | −0.362 | 0.097 | −3.751 | 0.001 | Treatment |
| Richness | ||||||||
| Treatment | Alpha.Score | A+ T- P- | A- T- P- | −0.064 | 0.046 | −1.388 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A- T- P- | −0.089 | 0.048 | −1.856 | 0.245 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T- P- | 0.067 | 0.057 | 1.174 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A- T+ P- | A+ T- P- | −0.025 | 0.050 | −0.502 | ≥0.25 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A+ T- P- | 0.131 | 0.059 | 2.226 | 0.115 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T+ P- | 0.156 | 0.061 | 2.586 | 0.047 | Treatment |
| Significant Results:Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE | ||||||||
| Tukey(Displacement ~ Treatment); All treatments at 60 DPE | ||||||||
| term | .y. | group1 | group2 | estimate | std.error | statistic | adj.p.value | Variable |
|---|---|---|---|---|---|---|---|---|
| Shannon | ||||||||
| Treatment | Alpha.Score | A+ T+ P- | A- T- P- | −0.191 | 0.067 | −2.860 | 0.022 | Treatment |
| Simpson | ||||||||
| Treatment | Alpha.Score | A+ T+ P- | A- T- P- | −0.332 | 0.092 | −3.627 | 0.002 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A+ T- P- | −0.355 | 0.094 | −3.765 | <0.001 | Treatment |
| Treatment | Alpha.Score | A+ T+ P- | A- T+ P- | −0.362 | 0.097 | −3.751 | 0.001 | Treatment |
| Richness | ||||||||
| Treatment | Alpha.Score | A+ T+ P- | A- T+ P- | 0.156 | 0.061 | 2.586 | 0.047 | Treatment |
# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressNoParaExp"
# GLM Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]] <-
test.beta.all_metrics$volatility %>%
# Filter
dplyr::filter(Treatment %in% c(
"A- T- P-",
# "A- T- P+",
"A+ T- P-",
# "A+ T- P+",
"A- T+ P-",
# "A- T+ P+",
"A+ T+ P-"#,
# "A+ T+ P+"
)) %>%
# Clean up cell value names by removing any strings including and after "__"
cutCellNames(col = "Beta.Metric", sep = "__") %>%
# Filter out raw data, leave only normalized data
dplyr::filter(Data.Type != "Raw") %>%
# Run GLM models
run_glm_models_beta(formula_str = "Displacement ~ Treatment",
beta_metric_col = "Beta.Metric",
beta_score_col = "Displacement",
family_str = "gaussian")
# ANOVA Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["ANOVA"]] <-
run_glm_anova_beta(beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]])
# Tukey Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["Tukey"]] <-
test.beta.all_metrics$volatility %>%
# Filter
dplyr::filter(Treatment %in% c(
"A- T- P-",
# "A- T- P+",
"A+ T- P-",
# "A+ T- P+",
"A- T+ P-",
# "A- T+ P+",
"A+ T+ P-"#,
# "A+ T+ P+"
)) %>%
# Clean up cell value names by removing any strings including and after "__"
cutCellNames(col = "Beta.Metric", sep = "__") %>%
# Filter out raw data, leave only normalized data
dplyr::filter(Data.Type != "Raw") %>%
# Run Tukey test on GLM models
run_tukey_glm_beta(., "Displacement", "Beta.Metric", c("Treatment"),
family_str = "gaussian",
group_by_var = NULL)
| GLM Results - Beta | |||||
| glm(Displacement ~ Treatment); Prior stressors and no parasite exposure | |||||
| term | estimate | std.error | statistic | p.value | p.adj.sig |
|---|---|---|---|---|---|
| Bray-Curtis_norm | |||||
| (Intercept) | −0.061 | 0.027 | −2.279 | 0.023 | * |
| A+ T- P- | 0.016 | 0.039 | 0.414 | ≥0.25 | ns |
| A- T+ P- | 0.081 | 0.039 | 2.059 | 0.040 | * |
| A+ T+ P- | −0.020 | 0.041 | −0.475 | ≥0.25 | ns |
| Canberra_norm | |||||
| (Intercept) | −0.176 | 0.023 | −7.610 | <0.001 | **** |
| A+ T- P- | 0.040 | 0.033 | 1.213 | 0.226 | ns |
| A- T+ P- | 0.179 | 0.034 | 5.285 | <0.001 | **** |
| A+ T+ P- | 0.155 | 0.035 | 4.359 | <0.001 | **** |
| ANOVA of GLM - Beta | ||||
| ANOVA(GLM(Displacement ~ Treatment), type = 2); Prior stressors and no parasite exposure | ||||
| term | statistic | df | p.value | sig |
|---|---|---|---|---|
| Bray-Curtis_norm | ||||
| Treatment | 6.694 | 3.000 | 0.082 | ns |
| Canberra_norm | ||||
| Treatment | 38.039 | 3.000 | <0.001 | **** |
| Pairwise Tukey's HSD, p.adj: Dunnett - Beta | ||||||||
| Tukey(Displacement ~ Treatment); Prior stressors and no parasite exposure | ||||||||
| term | .y. | group1 | group2 | estimate | std.error | statistic | adj.p.value | Variable |
|---|---|---|---|---|---|---|---|---|
| Bray-Curtis_norm | ||||||||
| Treatment | Beta.Score | A+ T- P- | A- T- P- | 0.016 | 0.039 | 0.414 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P- | A- T- P- | 0.081 | 0.039 | 2.059 | 0.166 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A- T- P- | −0.020 | 0.041 | −0.475 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P- | A+ T- P- | 0.065 | 0.040 | 1.621 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A+ T- P- | −0.036 | 0.042 | −0.851 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A- T+ P- | −0.101 | 0.042 | −2.371 | 0.083 | Treatment |
| Canberra_norm | ||||||||
| Treatment | Beta.Score | A+ T- P- | A- T- P- | 0.040 | 0.033 | 1.213 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P- | A- T- P- | 0.179 | 0.034 | 5.285 | <0.001 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A- T- P- | 0.155 | 0.035 | 4.359 | <0.001 | Treatment |
| Treatment | Beta.Score | A- T+ P- | A+ T- P- | 0.138 | 0.034 | 4.015 | <0.001 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A+ T- P- | 0.114 | 0.036 | 3.164 | 0.009 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A- T+ P- | −0.024 | 0.037 | −0.666 | ≥0.25 | Treatment |
# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressNoParaExp_60DPE"
# GLM Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]] <-
test.beta.all_metrics$displacement %>%
dplyr::filter(Treatment %in% c(
"A- T- P-",
# "A- T- P+",
"A+ T- P-",
# "A+ T- P+",
"A- T+ P-",
# "A- T+ P+",
"A+ T+ P-"#,
# "A+ T+ P+"
)) %>%
dplyr::filter(Time == 60) %>%
# Clean up cell value names by removing any strings including and after "__"
cutCellNames(col = "Beta.Metric", sep = "__") %>%
# Filter out raw data, leave only normalized data
dplyr::filter(Data.Type != "Raw") %>%
# Run GLM models
run_glm_models_beta(formula_str = "Displacement ~ Treatment",
beta_metric_col = "Beta.Metric",
beta_score_col = "Displacement",
family_str = "gaussian")
# ANOVA Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["ANOVA"]] <-
run_glm_anova_beta(beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]])
# Tukey Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["Tukey"]] <-
test.beta.all_metrics$displacement %>%
dplyr::filter(Treatment %in% c(
"A- T- P-",
# "A- T- P+",
"A+ T- P-",
# "A+ T- P+",
"A- T+ P-",
# "A- T+ P+",
"A+ T+ P-"#,
# "A+ T+ P+"
)) %>%
dplyr::filter(Time == 60) %>%
# Clean up cell value names by removing any strings including and after "__"
cutCellNames(col = "Beta.Metric", sep = "__") %>%
# Filter out raw data, leave only normalized data
dplyr::filter(Data.Type != "Raw") %>%
# Run Tukey test on GLM models
run_tukey_glm_beta(., "Displacement", "Beta.Metric", c("Treatment"),
family_str = "gaussian",
group_by_var = NULL)
| GLM Results - 60 DPE | |||||
| glm(Displacement ~ Treatment); Prior stressors and no parasite exposure at 60 DPE | |||||
| term | estimate | std.error | statistic | p.value | p.adj.sig |
|---|---|---|---|---|---|
| Bray-Curtis_norm | |||||
| (Intercept) | −0.060 | 0.043 | −1.398 | 0.165 | ns |
| A+ T- P- | 0.007 | 0.064 | 0.109 | ≥0.25 | ns |
| A- T+ P- | 0.043 | 0.067 | 0.640 | ≥0.25 | ns |
| A+ T+ P- | −0.119 | 0.079 | −1.501 | 0.136 | ns |
| Canberra_norm | |||||
| (Intercept) | −0.124 | 0.040 | −3.122 | 0.002 | ** |
| A+ T- P- | −0.089 | 0.059 | −1.498 | 0.137 | ns |
| A- T+ P- | 0.099 | 0.062 | 1.602 | 0.112 | ns |
| A+ T+ P- | 0.123 | 0.073 | 1.680 | 0.096 | ns |
| ANOVA of GLM - 60 DPE | ||||
| ANOVA(GLM(Displacement ~ Treatment), type = 2); Prior stressors and no parasite exposure at 60 DPE | ||||
| term | statistic | df | p.value | sig |
|---|---|---|---|---|
| Bray-Curtis_norm | ||||
| Treatment | 3.890 | 3.000 | ≥0.25 | ns |
| Canberra_norm | ||||
| Treatment | 12.030 | 3.000 | 0.007 | ** |
| Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE | ||||||||
| Tukey(Displacement ~ Treatment); Prior stressors and no parasite exposure at 60 DPE | ||||||||
| term | .y. | group1 | group2 | estimate | std.error | statistic | adj.p.value | Variable |
|---|---|---|---|---|---|---|---|---|
| Bray-Curtis_norm | ||||||||
| Treatment | Beta.Score | A+ T- P- | A- T- P- | 0.007 | 0.064 | 0.109 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P- | A- T- P- | 0.043 | 0.067 | 0.640 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A- T- P- | −0.119 | 0.079 | −1.501 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P- | A+ T- P- | 0.036 | 0.069 | 0.514 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A+ T- P- | −0.126 | 0.082 | −1.545 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A- T+ P- | −0.162 | 0.084 | −1.931 | 0.213 | Treatment |
| Canberra_norm | ||||||||
| Treatment | Beta.Score | A+ T- P- | A- T- P- | −0.089 | 0.059 | −1.498 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P- | A- T- P- | 0.099 | 0.062 | 1.602 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A- T- P- | 0.123 | 0.073 | 1.680 | ≥0.25 | Treatment |
| Treatment | Beta.Score | A- T+ P- | A+ T- P- | 0.187 | 0.064 | 2.921 | 0.018 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A+ T- P- | 0.212 | 0.076 | 2.804 | 0.025 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A- T+ P- | 0.025 | 0.078 | 0.321 | ≥0.25 | Treatment |
| Significant Results:Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE | ||||||||
| Tukey(Displacement ~ Treatment); All treatments at 60 DPE | ||||||||
| term | .y. | group1 | group2 | estimate | std.error | statistic | adj.p.value | Variable |
|---|---|---|---|---|---|---|---|---|
| Canberra_norm | ||||||||
| Treatment | Beta.Score | A- T+ P- | A+ T- P- | 0.187 | 0.064 | 2.921 | 0.018 | Treatment |
| Treatment | Beta.Score | A+ T+ P- | A+ T- P- | 0.212 | 0.076 | 2.804 | 0.025 | Treatment |
# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressNoParaExp_60DPE"
#### Dispersion --------------------------------------------------------------
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]] <-
run_BetaDispersion(dist.matrix = beta.dist.mat[[tmp.resSubSection]],
var = c("treatment_group"))
##### ANOVA ------------------------------------------------------------------
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA"]] <-
get_HoD_anova(betaDisper = beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]],
var = c("treatment_group"))
### Table
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA.Table"]] <-
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA"]] %>%
# Create the table
dplyr::group_by(Beta.Metric) %>%
set_GT(var = "p.value", group.by = "Beta.Metric") %>%
# Title/caption
gt::tab_header(
title = "ANOVA: Homogeneity of Dispersion",
subtitle = "ANOVA(Beta Disperson ~ Treatment); Prior stressors and no parasite exposure at 60 DPE"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "lightgreen"),
locations = gt::cells_body(
columns = c("p.value", "sig"),
rows = p.value < 0.05
)
)
##### Tukey ------------------------------------------------------------------
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey"]] <-
get_HoD_tukey(betaDisper = beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]],
var = c("treatment_group"))
### Table
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey.Table"]] <-
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey"]] %>%
dplyr::mutate(
group1 = case_when(
group1 == "Control" ~ "A- T- P-",
group1 == "Parasite" ~ "A- T- P+",
group1 == "Antibiotics" ~ "A+ T- P-",
group1 == "Antibiotics_Parasite" ~ "A+ T- P+",
group1 == "Temperature" ~ "A- T+ P-",
group1 == "Temperature_Parasite" ~ "A- T+ P+",
group1 == "Antibiotics_Temperature" ~ "A+ T+ P-",
group1 == "Antibiotics_Temperature_Parasite" ~ "A+ T+ P+",
),
group2 = case_when(
group2 == "Control" ~ "A- T- P-",
group2 == "Parasite" ~ "A- T- P+",
group2 == "Antibiotics" ~ "A+ T- P-",
group2 == "Antibiotics_Parasite" ~ "A+ T- P+",
group2 == "Temperature" ~ "A- T+ P-",
group2 == "Temperature_Parasite" ~ "A- T+ P+",
group2 == "Antibiotics_Temperature" ~ "A+ T+ P-",
group2 == "Antibiotics_Temperature_Parasite" ~ "A+ T+ P+",
)) %>%
# Create the table
dplyr::group_by(Beta.Metric) %>%
set_GT(var = "adj.p.value", group.by = "Beta.Metric") %>%
# Title/caption
gt::tab_header(
title = "Tukey: Homogeneity of Dispersion",
subtitle = "Tukey(Beta Disperson ~ Treatment); Prior stressors and no parasite exposure at 60 DPE"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "lightgreen"),
locations = gt::cells_body(
columns = c("adj.p.value", "sig"),
rows = adj.p.value < 0.05
)
)
| ANOVA: Homogeneity of Dispersion | ||||||
| ANOVA(Beta Disperson ~ Treatment); Prior stressors and no parasite exposure at 60 DPE | ||||||
| term | df | sumsq | meansq | statistic | p.value | sig |
|---|---|---|---|---|---|---|
| bray | ||||||
| treatment_group | 3.000 | 0.289 | 0.096 | 3.589 | 0.016 | * |
| Residual | 103.000 | 2.767 | 0.027 | NA | NA | NA |
| canberra | ||||||
| treatment_group | 3.000 | 0.076 | 0.025 | 6.573 | <0.001 | *** |
| Residual | 103.000 | 0.395 | 0.004 | NA | NA | NA |
| Tukey: Homogeneity of Dispersion | ||||||||
| Tukey(Beta Disperson ~ Treatment); Prior stressors and no parasite exposure at 60 DPE | ||||||||
| .y. | term | group1 | group2 | estimate | conf.low | conf.high | adj.p.value | sig |
|---|---|---|---|---|---|---|---|---|
| bray | ||||||||
| Distance | treatment_group | A+ T- P- | A- T- P- | −0.036 | −0.142 | 0.070 | ≥0.25 | ns |
| Distance | treatment_group | A- T+ P- | A- T- P- | −0.064 | −0.174 | 0.046 | ≥0.25 | ns |
| Distance | treatment_group | A+ T+ P- | A- T- P- | −0.162 | −0.294 | −0.031 | 0.009 | ** |
| Distance | treatment_group | A- T+ P- | A+ T- P- | −0.028 | −0.143 | 0.087 | ≥0.25 | ns |
| Distance | treatment_group | A+ T+ P- | A+ T- P- | −0.126 | −0.262 | 0.009 | 0.077 | ns |
| Distance | treatment_group | A+ T+ P- | A- T+ P- | −0.098 | −0.237 | 0.040 | ≥0.25 | ns |
| canberra | ||||||||
| Distance | treatment_group | A+ T- P- | A- T- P- | −0.056 | −0.096 | −0.016 | 0.002 | ** |
| Distance | treatment_group | A- T+ P- | A- T- P- | 0.010 | −0.032 | 0.051 | ≥0.25 | ns |
| Distance | treatment_group | A+ T+ P- | A- T- P- | −0.023 | −0.072 | 0.027 | ≥0.25 | ns |
| Distance | treatment_group | A- T+ P- | A+ T- P- | 0.066 | 0.023 | 0.109 | <0.001 | *** |
| Distance | treatment_group | A+ T+ P- | A+ T- P- | 0.034 | −0.017 | 0.085 | ≥0.25 | ns |
| Distance | treatment_group | A+ T+ P- | A- T+ P- | −0.032 | −0.085 | 0.020 | ≥0.25 | ns |
# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressNoParaExp_60DPE"
#### Capscale ----------------------------------------------------------------
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.mod"]] <-
run_capscale(ps.cleaned_alpha %>% # is the phyloseq object we are looping over
microViz::ps_filter(Treatment %in% c(
"A- T- P-",
# "A- T- P+",
"A+ T- P-",
# "A+ T- P+",
"A- T+ P-",
# "A- T+ P+",
"A+ T+ P-"#,
# "A+ T+ P+"
)) %>%
microViz::ps_filter(Time == 60),
dist.matrix = beta.dist.mat[[tmp.resSubSection]],
formula_str = "dist ~ Treatment")
##### ADONIS ------------------------------------------------------------------
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS"]] <-
run_cap_adonis(ps.cleaned_alpha %>% # is the phyloseq object we are looping over
microViz::ps_filter(Treatment %in% c(
"A- T- P-",
# "A- T- P+",
"A+ T- P-",
# "A+ T- P+",
"A- T+ P-",
# "A- T+ P+",
"A+ T+ P-"#,
# "A+ T+ P+"
)) %>%
microViz::ps_filter(Time == 60),
dist.matrix = beta.dist.mat[[tmp.resSubSection]],
formula_str = "dist ~ Treatment",
by.method = "terms")
## Warning in tidy.anova(.): The following column names in ANOVA output were not
## recognized or transformed: SumOfSqs, R2
## Warning in tidy.anova(.): The following column names in ANOVA output were not
## recognized or transformed: SumOfSqs, R2
### Table
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS.Table"]] <-
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS"]] %>%
set_GT(var = "p.value", group.by = "Beta.Metric") %>%
# Title/caption
gt::tab_header(
title = "ADONIS2",
subtitle = "adonis2(Beta Distance ~ Treatment); Prior stressors and no parasite exposure at 60 DPE"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "lightgreen"),
locations = gt::cells_body(
columns = c("p.value", "sig"),
rows = p.value < 0.05
)
)
| ADONIS2 | ||||||
| adonis2(Beta Distance ~ Treatment); Prior stressors and no parasite exposure at 60 DPE | ||||||
| term | df | SumOfSqs | R2 | statistic | p.value | sig |
|---|---|---|---|---|---|---|
| bray | ||||||
| Treatment | 3.000 | 1.688 | 0.102 | 3.887 | 0.001 | *** |
| Residual | 103.000 | 14.910 | 0.898 | NA | NA | NA |
| Total | 106.000 | 16.598 | 1.000 | NA | NA | NA |
| canberra | ||||||
| Treatment | 3.000 | 3.157 | 0.086 | 3.234 | 0.001 | *** |
| Residual | 103.000 | 33.522 | 0.914 | NA | NA | NA |
| Total | 106.000 | 36.679 | 1.000 | NA | NA | NA |